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1.0  ABSTRACT 


This  report  describes  an  attempt  to  design  by  analytic  means  a class  of 
axlsymmetrlc  bodies  having  low  drag  in  incompressible  flow  for  the  case  where 
the  boundary  layer  is  fully  turbulent  over  the  entire  body.  The  drag  is  to 
be  made  small  solely  by  proper  shaping  of  the  body,  and  the  drag  coefficient  to 
be  ..i1n1m1*eu  is  that  based  on  a reference  area  equal  to  the  two-thirds  power 
of  the  body's  volume.  A comparison  of  various  drag-calculation  methods  shows 
that  the  Truckenbrodt  formula,  which  expresses  drag  as  an  Integral  of  a power 
of  the  potential -flow  surface  velocity,  is  sufficiently  accurate,  and  this 
formula  is  adopted  as  the  chief  analytic  tool.  Among  the  consequences  of  this 
choice  is  that  bodies  with  low  drag  at  one  Reynolds  number  have  low  drag  at 
all.  Drag  performances  of  various  types  of  bodies  are  compared,  including 
those  of  "cavitation  shapes"  derived  from  a new  inverse  potential -flow  program. 
Two-dimensional  optimization  studies  are  carried  out  using  an  integral  formula 
for  drag  analagous  to  that  of  Truckenbrodt.  The  principal  conclusion  is  that 
the  drag  coefficient  is  not  very  sensitive  to  body  shape  and  thus  that  no 
significant  drag  reductions  can  be  obtained  from  shaping  alone. 
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4.0  PRINCIPAL  NOTATION 


area  enclosed  by  the  profile  curve  of  a two-dimensional  body 

drag  coefficient;  based  on  two-thirds  power  of  volume  for  axlsymmetrlc 
bodies;  based  on  square  root  of  A for  two-dimensional  bodies 

skin  friction  coefficient 

drag;  used  In  some  figures  for  body  diameter 

ratio  of  boundary-layer  displacement  thickness  to  momentum  thickness 
Integral  defined  by  (13) 

Integral  defined  by  (19) 
body  length 

Reynolds  number  based  on  body  length 
radial  coordinate,  distance  from  symmetry  axis 
reference  area  equal  to  two-thirds  power  of  body  volume 
potential-flow  velocity  on  the  body  surface 
free-stream  velocity 

velocity  component  in  boundary  layer.  On  the  body  it  is  parallel  to 
the  local  body  surface.  In  the  wake  It  Is  parallel  to  the  free-stream 
direction 

coordinate  parallel  to  free-stream  direction.  For  axlsymmetrlc  flow 
the  x-axis  Is  the  symmetry  axis. 

In  section  6.0;  denotes  distance  in  the  boundary  layer  measured  normal 
to  the  body  surface.  In  section  7.1  and  In  all  two-dimensional  analysi 
It  denotes  a coordinate  perpendicular  to  the  free-stream  direction 

slope  angle  of  the  profile  curve  of  the  body  with  respect  to  the  x-axis 

boundary-layer  thickness 

boundary-layer  displacement  thickness 

boundary-layer  momentum  thickness 

fluid  density  - a constant 

momentum  area 

momentum  area  of  the  wake  far  downstream 
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5.0  INTRODUCTION  AND  SUMMARY 


This  report  presents  results  of  a study  whose  object  has  been  the 
design  of  axl symmetric  bodies  having  low  drag  at  high  Reynolds  nunber.  More 
precisely  the  problem  may  be  stated  as  follows.  Consider  an  axisymmetrlc  body 
at  zero  angle  of  attack  In  an  Incompressible  flow  whose  Reynolds  number  Is  so 
large  that  the  boundary  layer  on  the  body  may  b«  taken  fully  turbulent  over  the 
entire  length  of  the  body.  Let  the  drag  coefficient  be  based  on  a reference  area 
defined  as  the  two-thirds  power  of  the  body's  volume.  It  Is  desired  to  mini- 
mize this  drag  coefficient  solely  by  proper  shaping  of  the  body  without  resort 
to  suction,  compliant  surfaces,  additives,  etc.  Moreover,  this  minimization 
should  be  performed  for  a range  of  Reynolds  number.  While  the  problem  is  easy 
to  state,  it  Is  difficult  to  formulate  in  a manner  that  admits  of  fruitful 
analysis.  Accordingly,  the  present  work  has  consisted  of  several  Investiga- 
tions designed  to  shed  light  on  the  mechanisms  leading  to  low  drag  and  also 
to  determine  a number  of  bodies  having  as  low  a drag  as  possible.  The  chief 
emphasis  has  been  on  the  axisymmetric  problem,  but  the  related  two-dimensional 
problem  has  been  considered  also.  The  latter  Is  more  tractable  in  some  ways, 
and  the  results  obtained  have  definite  Implications  for  the  axisymmetric  case. 

The  first  part  of  the  study  concerns  methods  of  computing  the  drag  of  an 
axisymmetrlc  body  (section  6.0).  The  ideas  and  formulas  for  calculating  drag 
from  a momentum  deficit  approach  are  presented.  These  formula*-  are  Imple- 
mented with  two  different  types  of  boundary- layer  calculation  methods  and  also 
with  a set  of  simplifying  issumptlons  that  eliminates  the  need  for  an  explicit 
boundary-layer  calculation  and  calculates  drag  from  an  integral  over  the  body 
of  a power  of  the  potential  flow  velocity.  Drags  computed  by  all  three  mpthods 
are  compared  with  experimental  drags  obtained  for  a set  of  bodies  In  water. 

It  is  concluded  that  the  simplified  Integral  method  of  computing  drag  is 
sufficiently  accurate  and  so  much  simpler  than  the  others  that  it  is  used  In 
all  subsequent  calculations.  A particularly  desirable  feature  of  this  method 
Is  the  fact  that  the  Reynolds  number  enters  only  in  a multiplicative  factor. 
Thus,  if  a particular  body  has  the  lowest  drag  at  one  Reynolds  number,  it  has 
the  lowest  drag  at  all  Reynolds  numbers,  and  the  investigation  need  not  concern 
itself  with  Reynolds  number  directly. 
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The  drag  performance  of  a wide  variety  of  bodies  has  been  compared  using 
the  simplified  Integral  drag  formula  (section  7.0).  The  drag  coefficient  for 
a general  prolate  spheroid  can  be  written  down  analytically.  The  curve  of 
drag  coefficient  versus  fineness  ratio  for  prolate  spheroids  has  proven  a very 
useful  device  for  correlating  drags  for  all  bodies,  because  all  "good"  shapes 
have  drags  that  plot  near  this  curve.  The  curve  has  a shallow  minimum,  and 
thus  there  Is  a preferred  range  of  fineness  ratio  from  about  three  to  four. 
Oddly,  the  simplified  integral  drag  formula  does  not  strictly  apply  to  prolate 
spheroids,  from  which  the  flow  is  presumed  to  separate  at  their  blunt  aft  ends. 
However,  the  separation  can  be  removed  by  a small  amount  of  aft  streamlining  or 
"boattailing"  with  a very  small  effect  on  drag.  In  any  case  the  unseparated 
"good"  bodies  have  relatively  low  drags  that  essentially  follow  the  curve. 

Unlike  certain  two-dimensional  applications,  there  is  no  information  on 
what  properties  of  the  body  shape,  its  pressure  distribution,  or  its  skin- 
friction  distribution  lead  to  low  drag.  Thus  the  use  of  inverse  boundary- layer 
or  potential-flow  methods  is  of  limited  value  in  the  present  problem.  Never- 
theless, some  investigation  along  this  line  has  been  carried  out.  In  partic- 
ular "cavitation"  bodies  having  extensive  regions  of  constant  pressure  were 
studied  on  the  ground  that  they  should  have  relatively  low  maximum  velocities 
and  thus  possibly  low  values  of  the  drag  integral.  Unfortunately,  this  proved 
not  to  be  true  (section  7.3). 

In  implementing  the  above  investigation  a new  method  of  solution  for  the 
axisymmetrlc  potential  flow  problem  was  discovered  and  developed.  This  solu- 
tion, which  is  applicaole  to  completely  general  shapes.  Is  based  on  a conformal 
transformation  of  the  body  to  a circle  and  obtains  the  solution  as  a series  of 
Chobyshev  and  Legendre  polynomials.  Both  direct  and  inverse  problems  can  be 
treated.  While  the  application  of  this  procedure  to  the  low-drag  problem 
appears  to  be  of  limited  usefulness,  the  development  of  such  a method  Is  felt 
to  be  important  In  Itself. 

Two-dimensional  analysis  has  been  used  In  section  8.0  to  attempt  to 
determine  shapes  that  minimize  the  drag  coefficient  calculated  from  the 
simplified  drag  integral.  For  this  study  the  two-dimensional  drag  coefficient 
has  been  based  on  the  square  root  of  the  area  enclosed  by  the  body  profile  as 
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the  closest  analogy  to  the  normalization  based  on  the  two-thirds  power  of  the 
volume  that  Is  used  In  the  axlsymmetrlc  case.  Slender-body  theory  Is  shown 
to  be  Inadequate.  General  solutions  based  on  conformal  transformations  yield 
bodies  having  minimum  drag  coefficients.  However,  the  drags  thus  obtained  are 
only  very  slightly  less  than  those  for  simple  bodies  selected  at  random. 
Probably  the  same  Is  true  for  the  axlsymmetrlc  problem. 
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6.0  DIRECT  DRAG  CALCULATION  FOR  AXISYMMETRIC  BODIES 


6.1  Theory  for  Drag  Calculation  Based  on  Momentum  Deficit 

Because  the  straightforward  calculation  of  drag  by  integration  of  surface 
pressure  and  skin  friction  turns  out  to  be  inaccurate  [1],  the  drag  of  a body 
must  be  obtained  by  considering  the  deficit  of  momentum  far  downstream  in  the 
wake.  The  basic  analysis,  which  is  given  in  [2],  is  merely  outlined  here. 

Consider  an  axisymmetric  body  at  zero  angle  of  attack  as  shown  in  Fig.  1. 
The  boundary  layer  on  the  body  continues  into  the  wake  and  far  downstream 
there  is  a deficit  of  momentum  due  to  viscous  retardation  of  the  flow  along 
the  body  surface.  From  momentum  considerations  the  drag  of  the  body  is  given 
by 
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D = 2ttP  J ufU^  — u)rdr  = 2ttp  J — u)rdr  (1) 

0 0 

where  p is  fluid  density  (a  constant  in  the  present  study),  is  free 
stream  velocity,  u is  the  local  velocity  component  parallel  to  the  symmetry 
axis  (x-axis),  r is  radial  distance  from  the  symmetry  axis,  and  6 is  the 
radius  of  the  wake.  The  integral  in  (1)  is  taken  in  an  x « constant  plane 
which  for  convenience  Is  taken  as  far  downstream.  The  momentum  area  of  the 
wake  far  downstream  is  defined  as 


so  that 
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the  drag  coefficient  is 
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where  S is  some  reference  area,  which  in  the  present  study  is  equal  to  the 
two-thirds  power  of  the  volume  of  the  body.  It  is  desired  to  compute  this 
drag  in  terms  of  quantities  evaluated  on  the  body,  and  this  is  where  the  first 
approximation  must  be  introduced. 


The  momentum  thickness  of  the  boundary  layer  at  any  point  of  the  body  is 

6 

9 = f S (]  ” ir)  1 1 + r cos  a)  dy 

0 

In  (5)  the  various  quantities  have  their  customary  definitions,  which  are 
somewhat  different  than  those  of  equation  (1),  which  applies  to  a location  in 
the  wake.  Specifically,  y denotes  distance  normal  to  the  body  surface,  <5 
is  boundary-layer  thickness,  u is  uhe  velocity  component  parallel  to  the 
surface,  and  U is  the  value  of  v for  y = 6.  Moreover,  r Is  radial 
distance  from  the  symmetry  axis  to  a point  on  the  body  surface,  and  a is  the 
local  slope  angle  of  the  surface  with  respect  to  the  symmetry  axis.  The 
momentum  area  of  the  boundary  layer  at  any  location  on  the  body  is  defined  as 

x = 2nre  (6) 


Finally  define 


where  6*  is  the  usual  boundary-1  aye*"  displacement  thickness 


(7) 
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5*  = J |l  - 7j-)|l  + £ cos  a | dy  (8) 
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Now  the  required  approximation  is  introducea.  Young's  analysis  [2]  leads  to  the 
relationship 


X~=  X( 


U {H+5>/2 

IT 

oo ' 


(9) 


where  x»  U,  and  H In  (9)  are  all  evaluated  at  a particular  location  on 
the  body  - usually  "near"  the  trailing  edge.  This  gives  for  the  drag 
coefficient 

C„  - (lf*s)n  < 


'D  ■ T 


Granville  [3]  has  an  alternate  form  of  (9)  which  gives  as  an  alternative  to 
(10)  the  form 

, . 4 irre  W X(V8)IH(17/8) 


However,  for  most  of  the  present  work  (10)  has  been  used.  The  various  quan- 
tities appearing  in  (10)  must  be  evaluated  by  means  of  a boundary- layer  calcu- 
lation. During  the  present  study  three  boundary- layer  methods  were  considered, 
one  of  which  evaluates  these  quantities  directly  in  terms  of  potential -flow 
quantities. 


6.2  Finite-Difference  Bcundary-Layer  Method 

The  most  general  and  elaborate  type  of  boundary- layer  method  is  that  which 
solves  a partial  differential  equation  for  flow  In  the  boundary  layer.  The 
numerical  solution  is  effected  by  some  type  of  finite-difference  scheme, 
usually  beginning  with  the  front  stagnation  point  and  proceeding  downstream 
along  tne  body.  The  calculation  terminates  if  separation  is  encountered  but 
othe*\/1se  traverses  the  entire  body.  In  the  present  study  of  low-drag  bodies 
it  is  presumed  that  all  good  bodies  will  indeed  have  unseparated  boundary 
layers,  so  that  the  above  Is  not  a serious  restriction.  In  fact  one  advantage 
of  the  finite-difference  method  is  that  it  predicts  separation  (or  Its  absence) 
which  some  other  methods  do  not.  For  a laminar  boundary  layer  such  a technique 
may  be  made  numerically  exact.  However,  for  a turbulent  boundary  layer  an 
empirical  or  semi -empirical  "closure"  condition  is  required  that  usually  takes 
the  form  of  an  "eddy-viscosity"  law.  This  assumption  renders  the  calculation 
approximate.  A variety  of  closure  conditions  have  been  proposed.  The  method 
used  in  this  study  Is  that  due  to  Cebeci  and  Smith.  Details  of  the  procedure 
are  given  in  [4]. 


According  to  the  theory  of  section  6.1,  the  drag  of  a body  Is  obtained 
by  evaluating  (10)  at  the  downstream  end  of  the  body.  However,  (10)  can  be 
evaluated  at  any  point  of  the  body.  If  this  Is  done  and  the  results  plotted 
versus  position  along  the  body,  certain  modifications  In  the  use  of  the  formula 
become  evident.  There  are  three  different  types  of  cases;  examples  of  which 
are  shown  In  figures  2,  3 and  4.  Figure  2 shows  the  unambiguous  case  where 
the  drag  coefficient  Is  a smooth  function  of  position  all  over  the  body,  and 
extrapolating  It  to  the  aft  end  of  the  body  presents  no  problem.  This  extrapo- 
lated value  Is  then  taken  to  be  the  drag  of  the  body.  More  common,  however, 
is  the  type  of  case  shown  in  Figure  3,  where  the  calculated  drag  coefficient 
rather  suddenly  begins  to  Increase  more  rapidly  a short  distance  ahc.td  of  the 
aft  end  of  the  body  — about  90%  of  body  length  In  figure  3.  Extrapolating 
the  calculated  curve  for  such  a case  to  the  aft  end  of  the  body  gives  a value 
of  drag  that  is  far  too  large.  The  most  effective  procedure  appears  to  con- 
sist of  extrapolating  to  the  aft  end  from  the  smooth  portion  of  the  curve 
upstream  of  the  "break"  as  shown  In  figure  3.  The  third  type  of  case  is  that 
shown  in  figure  4.  It  is  characterized  by  having  a calculated  drag  coefficient 
that  is  not  monotonlcally  increasing.  Instead  there  Is  a maximum  a short  distance 
ahead  of  the  aft  end  of  the  body,  and  the  drag  decreases  downstream  of  this 
point.  The  only  reasonable  procedure  Is  to  use  the  maximum  computed  value  of 
drag  as  the  drag  of  the  body. 

Despite  the  somewhat  striking  differences  in  the  procedure  for  estimating 
drag  in  the  three  cases,  there  Is  no  real  difference  In  the  accuracy  of  the 
results  as  judged  by  comparison  with  experimental  values.  The  "nice"  case  of 
figure  2 agrees  with  experiment  no  better  than  the  "difficult"  case  of  figure  3. 

6.3  Moment urn- Integral  Boundary-Layer  Method 

Momentum-lntegra1  boundary-layer  methods  are  characterized  by  the  fact 
that  they  use  certain  assumptions,  such  as  the  form  of  the  velocity  profile 
In  the  boundary  layer,  to  reduce  the  partial  differential  equation  of 
boundary-layer  flow  to  an  ordinary  differential  equation  (or  set  of  such 
equations)  along  the  body  surface.  Patel  [5]  has  developed  a method  of  this 
type  for  calculating  thick  axisymmetrlc  boundary  layers  and  Nakayama  and  Patel 
[6]  have  used  this  method  together  with  Granville's  formula  (11)  to  calcu- 
late drag.  As  in  the  third  type  of  case  of  section  6.2,  the  drag  calculated 
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in  this  way  does  not  increase  monotonlcally  along  the  body  hut  has  a maximum, 
which  is  taken  as  the  drag  of  the  body. 

6 . 4 Simplified  Integral  Drag  Calculation 

The  drag  calculation  methods  of  the  previous  two  sections  use  & boundary- 
layer  calculation  procedure  to  determine  the  quantities  that  enter  into  the 
drag  formula  (10).  The  boundary-layer  calculation  procedures  are  approximate 
but  are  still  fairly  elaborate  and  require  computer  programs  for  their  imple- 
mentation. Moreover,  the  drag  calculation  must  be  repeated  for  each  Reynolds 
number  of  interest.  By  making  certain  additional  approximations  the  quantities 
in  (10)  can  be  obtained  directly,  and  thus  the  need  for  a boundary- layer  calcu- 
lation method  can  be  eliminated  entirely.  In  the  present  study  this  has  been 
done  using  the  theory  of  Truckenbrodt  [7]. 


For  the  case  of  a fully  turbulent  boundary  layer  the  theory  of  [7]  gives 
the  momentum  thickness  at  a point  on  the  body  as 
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where  L is  body  length,  cf  is  total  skin-friction  coefficient,  and  where 
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(13) 


where  s denotes  arc  length  along  the  body  profile  curve,  and  where  the 
range  of  integration  is  from  the  front  stagnation  point  to  the  point  where 
9 is  to  be  evaluated.  If  this  is  used  in  (10)  the  drag  coefficient  is 
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Since  H > 1 It  Is  evident  that  If  the  potential-flow  velocity  is  used  in 
(14;,  then  Cg  = 0 at  the  aft  end  of  the  body.  However,  since  Cg  is  never 
negative.  It  obviously  has  a maxlmun  value,  which  generally  occurs  a short 
distance  ahead  of  the  aft  end  of  the  body  just  as  for  the  methods  of  sections 
6.2  and  6.3.  This  maxlmun  value  is  taken  as  the  drag  of  the  body. 
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Calculations  obtained  from  the  method  of  section  6.2  for  a large  number 
of  bodies,  have  shown  that  at  the  location  where  (14)  is  a maximun 

/ ij  v(H-l)/2 

0.975  < < 1.0  (15) 

for  all  bodies.  Thus  it  is  consistent  with  the  other  approximations  that  have 
been  made  to  set  this  quantity  equal  to  unity  and  to  write  the  drag  as 


(16) 


where  the  integral  is  over  the  entire  body  profile.  The  Reynolds  number  Re 
enters  only  in  the  skin-friction  coefficient  c^  which  the  Prandtl-Schlichting 
relation  [8]  gives  as 


0.455 

[logfRe)?755 


(17) 


Equation  (16)  Is  the  simplified  integral  drag  formula,  which  obtains 
drag  directly  from  the  potential -flow  velocity  distribution.  The  fact  that 
the  Reynolds  number  enters  only  in  a multiplicative  factor  means  that  if  one 
shape  has  lower  drag  than  another  at  one  Reynolds  number,  it  will  have  lower 
drag  at  all  Reynolds  numbers.  In  particular  an  optimum  shape, If  it  can  be 
found, is  Independent  of  Reynolds  number. 


6.5  Comparison  of  the  Methods 

The  relative  accuracies  of  the  three  methods  for  computing  drag  In  incompres 
sible  flow  have  been  evaluated  by  comparing  calculated  drags  with  the  highly 
accurate  experimental  results  obtained  by  Gertler  for  a series  of  bodies  of  revo- 
lution in  a towing  tank  [9].  For  these  comparisons  eight  of  the  20  bodies  of 
[9]  were  selected  in  such  a way  that  the  complete  ranges  of  fineness  ratio  and 
drag  values  are  well  represented. 

First  the  simplified  Integral  formula  (16)  is  compared  with  the  result 
of  using  the  finite- difference  boundary-layer  method  of  section  6.2  In  the 
basic  drag  formula  (10).  Figures  5 and  6 compare  these  two  methods  with 
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experiment  for  Reynolds  numbers  of  10  million  and  20  million,  respectively. 

The  simplified  Integral  formula  (16)  Is  definitely  more  accurate.  Its  RMS 
errors  are  approximately  60%  of  those  obtained  by  the  use  of  the  finite- 
difference  boundary-layer  method.  Formula  (16)  is  thus  much  to  be  preferred 
because  of  Its  simplicity. 

Next  the  momentum- integral  boundary-layer  method  of  section  6.3  used  with 
Granville’s  drag  formula  (11)  is  compared  with  the  simplified  Integral  formula 
(16)  and  with  the  result  of  using  the  finite-difference  boundary- layer  method 
of  section  6.2  in  the  basic  drag  formula  (10).  These  comparisons,  which  are 
shown  in  figures  7 and  8,  respectively,  are  for  the  same  set  of  bodies  from 
[9]  at  a Reynolds  number  of  20  million.  The  calculation  based  on  the  momentun- 
Integrai  is  definitely  the  mo!  accurate  of  the  three  methods  studied.  Its 
RMS  error  is  about  57*  of  that  obtained  by  the  simplified  integral  formula  and 
only  35%  of  that  obtained  using  the  finite-difference  boundary- layer  method. 

Despite  the  above  finding  the  simplified  integral  formula  (16)  has  been 
used  for  drag  calculation  in  this  study.  One  important  consideration  was 
that  the  momentum- integral  method  was  not  obtained  by  the  authors  until  rather 
late  in  the  study.  The  other  important  factor  is  the  greater  simplicity  of 
fotmula  (16),  including  the  fact  that  Reynolds  number  need  not  be  considered 
explicitly  in  the  investigation.  Despite  the  somewhat  inferior  accuracy  of 
formula  (16),  it  seemed  unlikely  that  any  overall  conclusions  drawn  on  the 
basis  of  this  formula  would  be  drastically  overturned  by  use  of  the  momentum- 
integral  method.  The  finite-difference  method  has  been  applied  after  the  fact 
to  the  best  bodies  to  verify  that  their  boundary  layers  are  unseparated. 


7.0  INVESTIGATION  OF  LOW-DRAG  SHAPES  USING  THE 
SIMPLIFIED  INTEGRAL  DRAG  CALCULATION 
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7 . 1 Comparison  of  Two-Dimensional  and  Axlsymmetrlc  Formulas 

Although  the  principal  Interest  In  the  present  study  Is  axlsymmetrlc 
flow,  two-dimensional  flow  has  been  considered  also,  because  It  Is  simpler 
from  some  standpoints  and  It  should  be  at  least  qualitatively  similar.  In 
fact  the  similarities  between  the  two-dimensional  and  axlsymmetrlc  cases  can 
most  easily  be  seen  by  examining  the  two  forms  of  the  simplified  drag  Integral. 
Spence  [10]  reduces  the  drag  of  a two-dimensional  body  to  an  Integral  of  the 
potential -flow  velocity  In  a manner  analogous  to  the  theory  of  Truckenbrodt, 
quoted  In  section  6.4,  for  axl symmetric  bodies.  The  two-dimensional  drag  coef- 
ficient must  be  normalized  with  respect  to  a reference  length  that  Is  as 
analagous  as  possible  to  the  reference  area  (two-thirds  power  of  volune)  used 
for  the  axlsymmetrlc  drag  coefficient.  The  length  chosen  Is  the  square  root 
of  the  area  enclosed  by  the  profile  curve  of  the  body.  With  these  choices 
the  drag  coefficient  for  both  axisynmetrlc  and  two-dimensional  bodies  may 
be  written 


(18) 


where  x and  y denote,  respectively,  distance  parallel  and  perpendicular  to 
the  free  stream.  As  before  U Is  surface  velocity,  s Is  arc  length,  and 
the  Integrals  are  over  the  entire  body.  The  exponents  are 
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7.2  Prolate  Spheroids  and  Pointed  Shapes 

The  function  to  be  adjusted  to  obtain  the  minimum  Is  of  course  the  body 
shape  y(x).  Some  direct  attacks  on  the  minimization  problem  are  outlined  In 
subsequent  sections.  In  this  section  attention  Is  concentrated  on  an  attempt 
to  find  good  shapes  simply  by  using  (16)  to  calculate  results  for  a variety  of 
bodies. 

The  drag  coefficient  of  (16)  can  be  obtained  analytically  for  prolate 
spheroids  as  a function  of  fineness  ratio.  This  function  Is  graphed  In 
Figure  9.  Drag  values  are  shown  for  a Reynolds  number  of  10  million,  but  at 
other  Reynolds  numbers  only  the  level  of  the  curve  Is  changed  not  Its  shape. 

It  can  be  seen  that  there  Is  a shallow  minimum  at  a fineness  ratio  somewhat 
greater  than  three.  This  curve  has  proved  quite  useful  In  analyzing  bodies, 
because  It  turns  out  that  on  the  basis  of  approximately  50  calculations  most 
good  bodies  give  drags  that  plot  very  near  this  curve.  It  might  be  argued 
that  since  prolate  spheroids  have  blunt  aft  ends  near  which  the  boundary 
layer  separates,  the  curve  of  figure  9 Is  not  meaningful.  To  Investigate  this 
several  prolate  spheroids  were  fitted  with  conical  boattalls.  Equation  (16) 
was  evaluated  for  these  bodies  and  the  results  lie  close  to  the  basic  curve 
as  shown  in  figure  9.  Absence  of  separation  was  verified  for  these  bodies 
by  finite-difference  boundary-layer  calculations.  By  way  of  example,  figure 
10  shows  results  for  a 40%-thlck  prolate  spheroid  that  was  fitted  with  a 
conical  boattall  to  produce  a body  of  fineness  ratio  3.25,  the  flow  about  which 
* Is  effectively  unseparated.  For  fineness  ratios  below  about  3,  absence  of 

separation  could  not  be  obtained,  and  it  Is  concluded  that  the  curve  of 
figure  9 Is  only  of  theoretical  Interest  to  the  left  of  this  point.  The  cal- 
culated drags  for  the  Gertler  bodies  (figure  5)  are  also  shown  In  figure  9 and 
they  also  lie  on  the  basic  curve.  The  calculations  Imply  a preferred  range 
of  fineness  ratio  from  3 to  4.  Experimental  results  may  alter  this  conclusion 
somewhat  — probably  In  the  direction  of  increased  fineness  ration. 

7.3  Constant-Pressure  Bodies  Derived  from  an  Inverse  Procedure 
I Since  the  Integrand  of  the  simplified  drag  Integral  (16)  depends  on  a 

higher  power  of  surface  velocity  than  of  local  body  ordinate.  It  was  hypothe- 
sized that  the  search  for  low  drag  bodies  might  better  be  conducted  In 
terms  of  velocity  distribution  rather  than  In  terms  of  shape  directly.  That 


is.  It  was  thought  that  certain  velocity  distributions  might  have  low  values 
of  drag  associated  with  the  corresponding  shapes.  However,  no  general  principles 
of  this  sort  could  be  formulated.  It  does  appear  from  the  form  of  (16)  that 
keeping  the  maximum  velocity  as  small  as  possible  should  be  desirable.  This 
leads  to  consideration  of  velocity  distributions  which  are  constant  at  their 
maximum  value  over  much  of  the  body,  i.e.,  "cavitation  shapes."  Unfortunately, 

It  turns  out  that  their  performance  is  not  much  different  from  prolate  spheroids 
(see  below). 

Suppose  a surface  velocity  distribution  can  be  selected  by  some  criterion. 
Determination  of  the  body  shape  corresponding  to  this  velocity  distribution 
requires  an  Inverse  potential -flow  method.  When  the  present  study  was  initiated 
the  development  of  such  a method  for  the  axi symmetric  case  was  an  unsolved 
problem.  Early  In  the  study  [11]  appeared,  which  solves  the  inverse  problem 
by  iterative  use  of  a direct  method.  As  part  of  the  present  stud.v  a new  type 
of  direct-and-inverse  solution  was  developed  for  the  problem  of  axisymmetric 
potential  flow.  The  method  of  solution  depends  on  conformal  transformation 
and  a series  solution  in  terms  of  Legendre  and  Chebyshev  polynomials.  A 
detailed  presentation  is  contained  in  Appendix  A. 

Based  on  the  reasoning  outlined  above,  "cavitation  shapes"  were  determined 
by  both  the  method  of  [11]  and  that  of  Appendix  A.  The  velocity  distributions 
that  were  called  for  had  extensive  regions  of  constant  velocity.  Because  of 
the  details  of  the  Inverse  procedures  the  bodies  that  resulted  had  velocity 
distributions  that  were  almost,  but  not  quite,  constant.  Figure  11  shows  a typ- 
• ical  result. 

Figure  12  shows  drag  coefficients  computed  by  (16)  for  nine  bodies  designed 
to  have  virtually  constant  maximum- velocity  regions.  Four  were  obtained  directly 
by  modifying  elliptic  contours.  The  other  five  are  products  of  the  inverse 
potential -flow  methods  mentioned  above.  Four  of  the  five  have  pointed  aft  ends 
similar  to  the  one  shown  In  Figure  11.  The  fifth,  which  has  the  lowest  drag  of 
all,  is  symmetric  fore-and-aft  and  has  a blunt  aft  end.  Flow  about  the  pointed 
shapes  Is  unseparated.  From  Figures  9 and  12  it  may  be  concluded  that  the  per- 
formance represented  by  the  curve  for  prolate  spheroids  can  actually  be  obtained 
but  that  shapes  with  significantly  lower  drag  are  difficult  to  find. 
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8.0  APPLICATION  OF  MINIMIZATION  TECHNIQUES  TO  THE  TWO-DIMENSIONAL 
LOW-DRAG  PROBLEM  BASED  ON  THE  SIMPLIFIED  INTEGRAL  FORMULA 

8.1  Statement  of  the  Problem 

One  of  the  principal  reasons  for  concentrating  01.  the  simplified  Integral 
formula  for  drag  Is  that  application  of  a rigorous  minimization  technique  to 
this  Integral  appears  to  be  at  least  a possibility.  According  to  this  formula 
a body  with  minimum  drag  at  one  Reynolds  number  has  minimum  drag  at  all  Reynolds 
numbers.  For  simplicity  the  Initial  work  on  minimization  vas  done  for  the  two- 
dimensional  problem.  The  results  obtained  were  not  sufficiently  encouraging  to 
justify  much  additional  effort  for  the  axlsymetrlc  case. 

The  problem  is  to  select  a shape  y(x)  In  such  a way  that  the  two- 

dimensional  form  of  (18)  Is  minimized.  Alternatively,,  the  function  y (x ) must 

be  selected  to  minimize  the  Integral 

J « D6/5  = (const)  j U4ds  (19) 

body 

subject  to  the  constraint  t? at  the  Integral 

A = j ydx  (20) 

body 

take  on  a specified  value.  This  is  not  a well-defined  calculus-of-variations 
problem.  The  difficulty  Is  that,  while  U depends  on  y In  the  sense  that 
given  a complete  body  shape  y(x)  then  U(x)  can  be  determined,  there  is  no 

relation  between  local  values  of  U and  y. 

8.2  Slender-Body  Theory 

For  sufficiently  slender  bodies  the  standard  aerodynamic  techniques  ylelo 
the  following  relationship 

a - / t , X/L  mi 

o Vi  - 0 - »)* 
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while  (19)  becomes 


1 

J « f U4dt  (22) 

0 

Now  standard  calculus-of-varlatlons  techniques  can  be  applied.  The  result  Is 


U = (const) 
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This  velocity  has  fore-and-aft  symmetry  with  a 1/6-power  singularity  at  the  ends. 
Clearly  the  theory  has  broken  down.  However,  to  investigate  the  possibility 
that  the  theory  has  Indicated  certain  desirable  features  for  low-drag,  calcula- 
tions were  performed  for  bodies  having  "saddleback"  velocity  distributions 
(fore-and-aft  peaks  of  finite  size  with  a lower  velocity  region  In  the  middle). 
These  bodies  consist  of  semicircles  at  the  front  and  rear  joined  by  a constant- 
thickness region,  whose  length  may  be  varied  to  give  different  fineness  ratios. 
Results  are  discussed  in  section  8.4  below  where  the  bodies  are  Identified  as 
"circle-flat"  bodies.  They  turned  out  to  have  relatively  high  drag.  Thus, 
slender-body  theory  was  abandoned. 

8.3  A General  Procedure  Based  on  Conformal  Mapping 

In  two-dimensional  potential  flow  the  complete  solution  for  both  direct 
and  inverse  problems  can  be  expressed  in  terms  of  the  coefficients  of  the  con- 
formal mapping  that  maps  the  body  In  question  to  the  unit  circle.  This  fact 
permits  a rather  general  "brute  force"  solution  to  the  minimization  problem. 


Let  z(c)  be  the  conformal  mapping  that  carries  a body  shape  in  the  z-plane 
to  the  unit  circle  In  the  c-plane.  Then  It  can  be  shown  that  the  mapping 
derivative  must  satisfy  a number  of  constraints  and  that  a general  expression 
consistent  with  these  constraints  can  be  written  down  quite  simply.  This  topic 
Is  discussed  In  more  detail  In  appendices  A and  B,  where  the  general  formula  In 
question,  viz. 
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Is  justified.  The  angles  t2  and  are  twice  the  cone  angles  of  the  fore  and 
aft  end  of  the  body,  respectively.  Equation  (24)  defines  the  mapping  coeffic- 
ients ap  which  must  be  real  for  application  to  bodies  of  revolution  by  reason 
of  the  symmetry.  Also  closure  of  the  profile  requires  that  a-j  = (t2  — t-j  )A. 


Now  both  the  surface  velocity  on  the  body  and  the  area  enclosed  by  the 
profile  can  be  expressed  in  terms  of  these  coefficients,  and  thus  so  can  the 
drag  coefficient  based  on  square  root  of  area  (18).  The  most  general  case  of 
(24)  has  not  been  Investigated  fully,  and  a detailed  account  of  the  various 
complications  is  the  subject  of  appendix  B.  However,  for  a blunt  body 
( t i * Tg  = *)  certain  Important  numerical  simplifications  occur  which  enable 
this  case  to  be  computed  rather  slpply.  Specifically  the  drag  coefficient  for 
such  a body  can  be  written 


where 


and 


CD  = (const) 


sin  cudu 
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H(w,an)  = [(1  + a2  cos  2u  + a3  cos  3w  + 
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(27) 


A very  simple  numerical  searching  scheme  has  been  developed  for  finding  minimum 
values  of  a function  of  a finite  number  of  real  variables  which  has  an  obvious 
application  to  the  problem  of  minimizing  the  given  by  (25)  and  (26)  (see 
appendix  B,  section  B6).  Notice  that  this  Is  a minimum  problem  in  n variables, 
not  a calculus  of  variations  problem.  If  only  a2  Is  retained,  the  body  shape 
is  an  ellipse  and  the  procedure  gives  the  analytically  correct  thickness  ratio 
for  minimum  drag  coefficient.  Using  a large  nunber  of  terms  gives  the  result 
that  the  optimum  shape  is  actually  pointed,  because  the  body  obtained  has  its 
Infinite  slope  only  In  a very  small  region  of  the  fore  and  aft  ends,  i.e., 
the  procedure  is  trying  to  give  a pointed  shape.  Figure  13  shows  the  optimum 
"blunt"  body,  which  appears  pointed.  Moreover,  the  coefficients  an  of  (26) 
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for  the  optimum  body  appear  to  be  asymptotically  equal  to  those  of  an  expansion 
for  a particular  fractional  power.  The  indications  are  that  the  true  optimum 
has  an  equation  of  the  form 


which  is  the  simplest  form  of  (24),  namely  the  one  with  all  an  = 0.  The  equal 
fore  and  aft  angles  t appear  to  be  approaching  approximately  tt/3,  i.e.,  30° 
cone  angle  (see  appendix  B,  section  B7).  The  near-optimum  body  of  figure  13 
achieves  a drag  reduction  over  the  optimum  ellipse  of  0.15%  — a very  interesting 
and  discouraging  result. 

It  appears  that  by  using  the  new  type  of  axisynmetric  solution  described 
in  appendix  A,  an  analogous  minimization  procedure  could  be  applied  to  the 
axlsymmetric  case.  However,  all  the  indications  are  that  only  very  small 
reductions  in  drag  could  be  obtained. 

8.4  Miscellaneous  Two-Dimensional  Results 

While  the  emphasis  of the  present  study  is  on  axi symmetric  bodies,  some 
two-dimensional  cases  have  been  investigated  also.  Most,  but  not  all,  of  these 
cases  are  connected  with  the  optimization  studies  of  the  preceeding  sections. 

By  analogy  with  the  axisymmetrlc  case,  all  two-dimensional  bodies  considered 
have  a right-and-left  symmetry  with  respect  to  the  flow  direction  and  are, 
of  course,  nonlifting. 

Figure  14  summarizes  the  drag  coefficients  obtained  for  various  two- 
dimensional  bodies  from  the  integral  formula  (18).  The  curve  of  figure  14 
represents  the  analytic  expression  to  which  (18)  reduces  for  the  case  of 
ellipses.  The  drag  of  the  "near-optimum"  body  of  figure  13  is  shown,  as  well 
as  drags  for  a set  of  "circle-flat"  bodies  that  have  "saddleback"  velocity 
distributions.  As  discussed  in  section  8.2  these  bodies  were  suggested  by 
results  from  slender-body  theory,  but  it  can  be  seen  that  their  drags  are 
significantly  higher  than  ellipses. 

Figure  14  also  shows  drags  for  a series  of  symmetric  struts  that  were 
designed  to  have  zero  skin  friction  over  their  aft  portions  [12].  Despite 
this  their  total  drag  coefficients,  as  calculated  from  (18),  are  higher  than 
those  of  ellipses. 
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9.0  CONCLUSIONS 


In  computing  the  drag  of  an  ax i symmetric  body  with  fully  turbulent 
boundary  layer  in  incompressible  flow  the  momentum- integral  method  of 
Nakayama  and  Patel  [6]  Is  considerably  more  accurate  than  either  the 
method  based  on  a finite-difference  boundary- layer  calculation  [4]  or 
the  simplified  integral  formula  based  on  the  analysis  of  Trurkenbrodt  [7]. 
Of  the  latter  two  methods  the  simplified  integral  formula  is  more  accurate 
than  the  current  finite-difference  method  and  much  easier  to  use. 

If  the  simplified  integral  formula  is  used  to  compute  drag,  the  following 
conclusions  can  be  drawn  for  axisymmetric  bodies  with  fully  turbulent 
boundary  layers  in  incompressible  flow: 

a.  A shape  having  the  lowest  drag,  a one  Reynolds  number  has  the  lowest 
drag  at  all  Reynolds  numbers. 

b.  Shapes  with  fineness  ratios  in  the  range  three  to  four  have  the 
lowest  drag  coefficients  based  on  the  two-thirds  power  of  volume. 

c.  Drag  coefficient  is  insensitive  to  shape  and  no  shape  has  been  found 
with  significantly  lower  drag  than  a boattailed  prolate  spheroid. 

A more  accurate  drag  calculation  might  modify  these  conclusions 
slightly  but  would  probably  not  drastically  revise  them. 

The  two-dimensional  analogue  of  the  simplified  integral  formula  for  drag 
can  be  rigorously  optimized  in  terms  of  mapping  coefficients  to  find  a 
"near  optimum"  shape  that  supposedly  has  the  lowest  possible.  Unfortunately 
the  shape  so  determined  has  a drag  coefficient  almost  undetectably  less 
than  bodies  selected  at  random. 

Accurate  and  very  rapid  solutions  to  both  the  direct  and  the  inverse 
problems  of  axisymmetric  potential  flow  can  be  obtained  using  the  expan- 
sion method  of  appendix  A. 
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Figure  4.  Drag  computed  by  the  finite-difference  boundary- layer  method  as  a function  of  location 
for  a body  with  a cusped  aft  end. 


EXPERIMENTAL  DRAG  , FROM  [9] 


Figure  5.  Comparison  of  drags  computed  by  the  finite-difference  boundary- 
layer  method  and  the  simplified  integral  method  with  experimental 
data  for  a series  of  eight  bodies  at  a Reynolds  number  of  10  million. 


EXPERIMENTAL  DRAG,  FROM  [9] 


Figure  6.  Comparison  of  drags  computed  by  the  finite-difference  boundary- 
layer  method  and  the  simplified  Integral  method  with  experimental 
data  for  a series  of  eight  bodies  at  a Reynolds  number  of  20  million. 
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EXPERIMENTAL  DRAG  , FROM  [9] 


Figure  7.  Comparison  of  drags  computed  by  the  momentum- Integral  boundary-layer 
method  and  the  simplified  Integral  method  with  experimental  data  for 
a series  of  eight  bodies  at  a Reynolds  number  of  20  million. 


9 


O FINITE  DIFFERENCE 
A MOMENTUM  INTEGRAL 


RMS  ERRORS 
O.OOI74 
0.00062 


CALCULATED 


0.016  0 016  0.020  0.022  0.024 

EXPERIMENTAL  DRAG,  FROM  [9] 


0.026 


Figure  8.  Comparison  of  drags  computed  by  the  momentum-integral  boundary-layer 
method  and  the  finite-difference  boundary- layer  method  with  experi- 
mental data  for  a series  of  eight  bodies  at  a Reynolds  number  of 
20  million. 
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Figure  9.  Calculated  drag  coefficients  versus  fineness  ratio  at  a Reynolds  number  of  10  m.llion. 
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Calculated  drag  coefficients  versus  fineness  ratio  for  constant-pressure 
at  a Reynolds  number  of  10  million. 


e 13.  A two-dimensional  near-optimum  "blunt"  body. 
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APPENDIX  A 

A GENERAL  ANALYTICAL  METHOD  FOR  AXISYMMETRIC 
INCOMPRESSIBLE  POTENTIAL  FLOW 

ABSTRACT 

A method  is  presented  for  calculating  the  flow  field  about  bodies  of 
revolution  in  incompressible  potential  flow  as  a sequence  of  elementary  analytic 
functions  (Fourier,  Chebyshev  and  Legendre).  For  the  cases  of  greatest  inter- 
est to  practical  engineering  (cusp  and  blunt  trailing  edges)  comparison  with  an 
existing  highly  accurate  numerical  method  shows  that  convergence  is  good.  Only 
ten  terms  are  required  to  give  adequate  accuracy  for  negligible  computer  usage 
on  bodies  with  much  more  character  than  usually  encountered.  The  theory  can  be 
used  for  design  (inverse)  operation  to  produce  body  shapes  associated  in  some 
least  squares  sense  with  a desired  input  velocity.  It  is  shown  that,  when  used 
this  way,  again  good  results  are  obtained  with  remarkably  few  terms  or  iterative 
cycles. 


Al.  INTRODUCTION 

The  literature  of  axl symmetric  incompressible  potential  flow  abounds  with 
numerical  methods  of  the  direct  type.  In  which  Integral  equations  arc  solved 
eventually  by  matrix  methods  following  some  form  of  discretization.  This  char- 
acteristic traditional  path  starting  from  von  Karman's  [Al]  original  numerical 
method  and  arriving  at  such  sophisticated  schemes  as  the  higher-order  Neumann 
method  of  Hess  [A2]  has  not  only  spawned  no  general  analytical  techniques,  but 
even  very  few  special  solutions.  In  this  context  "analytical"  will  be  taken  to 
mean  solutions  in  terms  of  known  functions,  or  a sequence  of  known  functions, such 
that  the  convergence  and  computation  of  the  coefficients  allows  practically 
useful  calculations  to  be  carried  out  to  arbitrary  order.  Thus  solutions  for 
which  only  two  or  three  terms  can  be  obtained  (e.g.,  asymptotic  series  or 
matched  local  expansions,  etc.)  are  not  considered  to  be  analytical  In  this 
sense,  ihis  appendix  deals  with  the  development  of  a general  analytic  method 
for  axi symmetric  flow. 

In  two  dimensions  there  are  analytical  methods  for  isolated  bodies  (air- 
foils), and  experience  with  them  has  shown  that  they  have  certain  clear 
advantages  over  the  purely  numerical  schemes  - as  compensation  in  a very  real 
sense  for  loss  of  generality.  For  instance  — they  are  usually  quicker  and 
more  accurate;  they  frequently  provide  detailed  qualitative  insight  into  far 
field  or  special  local  behavior  (near  sharp  trailing  edges  as  an  example),  and 
they  often  have  a much  more  Immediate  application  to  inverse  problems.  In  the 
case  of  Isolated  airfoils  this  inverse  capacity.  In  which  profiles  are  designed 
for  a prescribed  velocity  distribution,  has  proved  very  valuable. 

One  of  the  most  Important  distinctions  between  the  two  dimensional  and 
axlsynmetric  cases  from  a theoretical  point  of  view  Is  that  the  former  flows 
can  be  treated  entirely  by  complex  variable  methods  Including  conformal  map- 
ping; whereas,  although  the  meridian  profile  of  an  axisymmetrlc  shape  Involves 
only  two  variables,  the  flow  field  Itself  is  not  Laplaclan  in  the  potential. 

This  means,  at  best,  that  one  could  perhaps  treat  the  geometry  by  mapping 
methods;  but  then  a much  more  sophisticated  differential  equation  would  have 
to  be  solved.  It  seems  that  this  largely  explains  the  almost  total  lack  of 
analytical  methods  even  though  the  actual  equations  of  motion  appear  decep- 
tively simple.  However,  one  island  in  the  sea  of  numerical  methods  is  the  work 
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| or  Kaplan  [A3],  who  did  use  a mapping  for  the  geometry  and  then  devised  an 

Iterative  scheme  for  the  velocity  potential.  Unfortunately  Kaplan*s  scheme 
rapidly  becomes  very  complicated  so  that  It  Is  only  feasible  to  work  out  the 
first  three  terms  In  detail,  and  quite  Impossible  to  program  a computer  to 
calculate  higher-order  terms  automatically. 

Nevertheless,  the  underlying  principle  of  this  calculation  was  a good  one. 
By  using  recent  advances  In  appreciation  of  conformal  mapping  methods  (James 
[A4])  It  Is  shown  below  that  a very  wide  class  of  bodies  can  be  parameterized 
In  a very  simple  way,  and  that,  consequently,  analytical  solutions  can  be  found 
for  general  axlsymmetrlc  shapes  by  a slightly  unexpected  series  assumption  - 
for  which  an  Indefinite  number  of  terms  can  be  computed  automatically.  From 
the  examples  given  It  Is  clear  th?t  convergence  can  be  expected  to  be  very 
good  In  general  when  used  to  calculate  the  direct  flow  about  a given  body. 
However,  since  this  Investigation  was  originally  stimulated  by  recent  interest 
[A5]  in  the  possibility  of  a design  (inverse)  procedure  for  axisymmetrlc  bodies. 
It  Is  even  more  significant  that  the  convergence  was  also  found  to  be  excellent 
on  the  cases  tested  when  the  theory  was  used  in  the  inverse  mode. 

A2.  EQUATIONS  OF  MOTION 

For  an  axlsymmetrlc  body  lying  In  the  (x,y)-plane  such  that  the  free- 
stream  at  infinity  is  parallel  to  the  axis  of  symmetry  (x),  the  velocity 
potential  (<f>)  satisfies  the  well-known  differential  aquation 

My&M y(yf£)‘°  (A,) 

with  the  additional  conditions  that  <f>  ^ x In  the  far  field  and  that  the 
derivative  of  <p  normal  to  the  profile  vanishes  on  the  profile.  It  is  quite 
standard  to  seek  orthogonal  curvilinear  transformations  x = xUj.  e2)» 
y = y(£j*  in  order  to  either  simplify  (Al)  or  the  profile  representation 
(thus  making  the  normal  derivative  condition  more  tractable).  Such  a trans- 
formation leads  to 
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where  and  h2  are  the  usual  metrics,  viz. 


However,  If  the  transformation  Is  In  particular  a conformal  mapping  from  the 
Si  + 1e2  plane  to  the  x + iy  plana  then  the  Cauchy-Rlemann  equations  hold, 
that  Is 

= 3v_  ax  s ay 

H1  h2*  3c2  Hl 

and  therefore  h(  = h2  which  reduces  (2)  to 


3 

/y  i4_\  + JL 

3*1 

PV  n2 

Even  though  (A3)  is  formally  no  simpler  than(Al)  It  Is  possible  to  find  general 
mappings  which  greatly  simplify  the  surface  boundary  condition  by  relating  the 
given  profile  to  one  of  a number  of  canonical  forms  In  an  auxiliary  (c^, 
plane.  Specifically  the  class  of  mappings  which  map  airfoil-like  profiles 
Into  the  unit  circle  have  been  found  to  be  of  great  utility  in  two  dimensions 
and  so  are  considered  further  for  this  application. 

A3.  UNIT  CIRCLE  MAPPING 

According  to  Rlemann’s  fundamental  mapping  theorem  any  profile  whose 
contour  encloses  a singly-connected  domain  (In  the  usual  engineering  sense 
of  e.g.  Woods  [A6])  can  be  mapped  Into  the  unit  circle.  However,  for  Airfoils 
and  axlsymmetrlc  bodies  of  Interest  to  engineering,  the  types  of  mapping  of 
sufficient  generality  and  simplicity  to  be  useful  are  quite  constrained.  The 
role  of  these  constraints  has  been  discussed  by  James  0M1  and  need  not  be 
elaborated  In  detail  here.  Briefly,  the  desired  characteristics  are  that: 

(a)  the  mapping  should  reflect  the  existence  of  a finite  tralllng-edge  angle 
discontinuity  (it  — t),  (b)  the  mapping  derivative  must  -*-1  as  jz|  -*•  <», 

(c)  the  profile  must  be  closed,  and  (d)  the  mapping  derivative  must  have 
neither  zeros  nor  singularities  on*  o-"  outside  the  unit  circle. 

*In  the  general  theory  of  mapping  methods  certain  boundary  singularities  may 
arise,  but  they  are  not  of  Interest  In  this  application. 
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I Let  ; ■ c + In  be  the  unit  circle  plane  then  It  Is  easy  to  show  that  a mapping 

I which  conforms  to  (a)  - (d)  Is  of  the  type 


(A4) 


where  g(c)  Is  an  analytic  function  with  neither  zeros  nor  singularities  out- 
side C (the  unit  circle)  and  which  -+1  as  |c|  -+  ».  Many  different  forms  of 
g(c)  can  be  used  to  generate  profiles,  but  for  the  present  purpose  the  more 
general  representation  of  g(c)  as  an  expansion 

9(C)  • 1 * (1  ~T/,)  +4  + ...  (A5) 

c 

Is  appropriate.  That  the  first  coefficient  a-|  Is  (1  - t/tt)  Is  a consequence 
of  the  closure  condition  (c). 
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It  Is  possible  to  prove  In  general  that  If  the  origin  In  the  z-plane  Is 
located  at  the  trailing  edge  then,  as  a consequence  of  (A4)  and  (A5), 


z 


:MJ 


,2-t/tt 


1 + a^S-j 


(A6) 


where  S2  ...  are  polynomials  of  a particularly  simple  form.  The  nature 
of  these  polynomials  Is  not  Important  here  since  (A 6)  will  not  be  u*ed  directly; 
however.  It  does  show  that  very  awkward  fractional  powers  will  arise  in  the 
analysis  of  the  differential  equation  (A3)  through  the  y factor  which  Is  the 
Imaginary  part  of  z.  In  order  to  test  the  utility  of  this  theory  it  was  first 
Investigated  for  the  cases  which  do  not  give  rise  to  fractional  powers,  namely 


r = n (blunt)  c g(c)  = 1 + + --j  + ...  (A7 ) 

r * 0 (cusp)  55  (l  — g(c)  - (l  ~ -^Wl  + ~ J (A8) 


and  the  remainder  of  this  study  deals  only  with  these  forms.  Fortunately,  the 
exclusion  of  fractional  powers  Is  not  very  Important  for  most  engineering 
applications,  because  those  bodies  of  practical  Interest  which  are  not  truly 
blunt  (like  ellipsoids)  usually  have  small  values  of  t/tt  for  which  the  rear 
stagnation  region  (where  the  velocity  falls  very  sharply  to  zero)  Is  very 
localized.  Hence,  It  Is  of  little  Importance  since  the  boundary  layer  tends 
to  "smudge"  the  details  very  close  to  the  trailing  edge  and  the  question  of  the 
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"proper"  combined  flow  In  this  region  Is  an  unsolved  problem  well  beyond  the 
scope  of  this  work. 


V 


v 


If  the  forms  1n(A7)  and(A8)  are  reorganized  and  Integrated  both  yield  a 
series 


z 


C + 


B1  B2 


(A9) 


where  C Is  an  arbitrary  real  constant  (again  choosing  the  trailing  edge  as 
the  origin  of  z)  and,  for: 

t ~ tt  (blunt),  B-j  — B2  “ a^/2,  ••• 

(A10) 

i ~0  (cusp)p  B«j  ■ dg  1 * B2  ” (flg  “ agz/Z,  ••• 

Naturally  in  dealing  with  the  unit  circle  In  the  c-plane  it  is  expedient  to 
use  the  angular  variable  (u)  so  that 

c = re1w  , r = 1 represents  the  profile 

but  z Is  not,  of  course,  an  analytic  function  of  r + 1u>,  so  to  circumvent 

this  it  is  necessary  to  use  the  variable 

r,  = e°  where  a = A + 1w,  r = eX  (All) 


For  axisymmetrlc  bodies  the  reflection  symmetry  ensures  that  all  the  an  and 
Bn  coefficients  are  real  so  that  using  (AH)  in  (A9)  gives  the  general  form  of 
y as 


with  the  boundary  conditions 


■ Bje”x  sin  to 

+ B2e_2x  sin  2w  + ... 

(A12) 

equation 

/y  +|-j 

fy  0 

(A13) 

y 3 A f 3o)  ' 

V aw/ 

^ ex  cos  u 

as  a 

(A14) 

= 0 

when  A = 0 
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Although  this  Is  a properly  posed  problem  which  has  so  far  defied  analytic, 
solution,  the  reduction  presented  In  this  section  Is  by  no  means  exhaustive; 

It  Is  merely  the  simplest.  Thus  the  question  of  convergence  has  been  tacitly 
Ignored.  Experience  has  shown  that  this  breakdown  Is  adequate  for  most  pur- 
poses, but  It  Is  clear  that  when  t = 0,  dz/de  has  a higher  order  zero  at 
the  trailing  edge  than  when  t = it.  If,  then,  convergence  Is  adequate  for  the 
blunt  cases  when  using  (A9)or  (A12)(wh1ch  effectively  disguises  the  zero  factor) 
it  will  be  less  so  for  the  cusp  cases.  Furthermore, (A9)  shows  quite  clearly 
that  one  should  remove  a zero  factor  from  the  expression  for  z In  both  cases 
for  even  better  convergence.  However,  It  is  not  possible  to  consider  all 
these  ramifications  here.  As  will  be  seen  the  basic  method  summarized  In  (A12), 
(A13)  and  (A  14) is  remarkably  good  and  behaves  exactly  as  these  conjectures  would 
lead  one  to  expect. 

A4.  SOME  REMARKS  ON  DIFFERENT  ITERATIVE  SCHEMES 

One  classical  method  of  solution  is  to  seek  a sequence  of  approximations 
4>0  • <t>2  such  that 


4 * *o  + *1  + *2  + *** 

and,  by  substitution  1nto(A13),  separate  a series  of  differential  equations 
each  of  which  can  be  solved  for  the  Individual  orders.  The  number  of  differ- 
ent ways  In  which  this  kind  of  analysis  can  be  effected  Is  quite  remarkable, 
but  the  varying  forms  of  the  successive  <|>n  functions  and  the  manner  and  rate 
of  their  convergence  is  even  more  surprising.  Since  the  natural  and  obvious 
assumptions  do  not  work.  It  would  be  remiss  ntt  to  Include  some  brief  discus- 
sion of  the  general  character  of  such  schemes. 


(a)  For  Instance, (A1 3)  can  be  expressed  In  the  form 


or  symbolically  as 


yv2*  + + y„»u  = 0 

yD$  + dyd4>  = 0 


(A15) 

(A16) 


Then  by  Identifying  the  successive  terms  of  y = y + y.  + ...  as 

\ _1L  0 1 

y = e sin  w,  y-j  * e sin  u>,  yg  s e sin  2u>  ...,  and  expanding  (A16) 
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one  Is  led  to  an  array  which  can  be  arranged  as: 


Vo  * Vi  +y0°2  + + dV*0 + dV*i + dy0d‘2  * ••• 

y10Q  + y^D2  + •••  + dy^d$^  ♦ dy^d$^  + dy^d$2  ^ •••  * 0 

y2D0  + y2Di  + y2°2  + • • • + dy2d*0  + dy2d*i  + dy2d*2  + • • • 

...  ...  (A17) 


The  simplest  scheme  which  successively  annihilates  the  left  hand  side  of  (A17) 
appears  to  be 


Do  “ 0 

Vi  ■ -dv*o 

V2  ■ ‘Vl  - dV*1  - dV4o 


(A18) 


and  It  Is  worth  noting  that  each  equation  Is  of  Poisson  form  and  so  readily 
solvable.  Furthermore,  although  (A18)  Is  not  unique  In  any  sense  because  any 
number  of  known  yn*s  could  be  appended  to  the  right  sides.  It  Is  at  least 
consistent  Inasmuch  as  each  order  does  not  Introduce  higher  harmonics  than 
required  for  the  Laplaclan  part.  That  1s,(A18)  could  well  be  solved  by  the 
assumptions 

<>o  = fQ(*)  cos  at,  $1  * V*)  cos  u,  $2  * cos  2(0»  ••• 

and,  proceeding  In  this  manner  with  the  proper  boundary  conditions,  readily 
leads  to 

4>  = 2 cosh  x cos  uj 
0 

» -(1  + x)e*A  cos  <i» 

$2  = |(1  + x + ^ x2)e  A — B-j  [(1  + 4x)e  A + e ^Aj|  cos  u» 

and  rapidly  Increasing  complexity.  It  Is  clear  that  this  scheme  suffers  from 
obvious  defects,  among  which  the  crucial  ones  are:  (1)  the  analysis  Increases 
rapidly  In  complexity  and  still  requires  solution  of  nonhomogeneous  differ- 
ential equations  at  each  step  and  so  cannot  be  mechanized,  (2)  the  convergence 
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In  a qualitative  sense  Is  poor  since  ^ has  only  just  Introduced  a B 
coefficient  - merely  the  first  - and  varies  only  as  cos  u>.  Part  of  the 
Inadequacy  of  this  scheme  arises  from  the  fact  that  each  equation  of  (A18) 
has  an  essentially  two-dimensional  character  so  that  the  "rate  of  Introduction 
of  axl symmetry"  Is  poor. 

(b)  Kaplan's  scheme  Is  much  more  sophisticated  and  has  an  axlsymmetrlc 
character  In  each  approximation  right  from  the  beginning.  It  can  be  represented 
In  terms  of  the  symbology  of  (A17)  by  the  sequence 

<*0  ♦ VDo  + (dyo  + dyl>d*o  * 0 

<yo  ♦ y,)D,  ♦ {dy0  ♦ dy,)d*,  - -y2D0  - dy2d*0  ^ 

(y0  ♦ y,  )02  + (dy0  + dy,  )dd2  - -y2D,  - y3D0  - dy2d*, 

-dy3d*o 

which  certainly  penetrates  the  array  faster  than  (A18).  In  fact  this  scheme 
not  only  Introduces  B2  and  cos  2u  In  «j>2»  but  the  first  term  Is  repre- 
sentative of  flow  around  an  ellipsoid,  not  merely  an  ellipse,  and  all  the  left- 
side operators  have  a more  nearly  axlsymmetrlc  character. 

Unfortunately  this  Improvement  In  the  character  of  each  term  has  been 
bought  at  the  expense  of  enormous  Increase  of  complexity.  Using  separation 
of  variables  on  each  equation  of  (A19)  leads  to  left  sides  which  give  rise  to 
two  Legendre  equations  Involving  the  Legendre  functions  of  both  kinds.  The 
function  of  second  kind  (Qn ) Is  awkward  to  handle, and  the  occurrence  of 
such  combinations  In  the  nonhomogeneous  part  of  each  equation  leads  to 
appalling  complexity  even  for  the  third  function  $3.  Thus  further  terms  by 
analysis  would  be  effectively  Impossible  and  computer  mechanization  unthink- 
able. However,  It  Is  Interesting  to  note  that,  even  so,  Kaplan's  approximation 
gave  a remarkably  good  answer  for  a Joukowskl  profile  of  revolution.  This 
provided  an  important  stimulus  to  the  work  reported  here,  and  Is  discussed 
further  in  paragraph  8 c . Kaplan's  method  in  general  also  made  it  clear 
that  1 teratl ve/analytic  methods  of  the  type  discussed  above  could  never  lead 
to  general  analytical  solutions  unless  some  means  could  be  found  to  make  each 
succeeding  order  depend  only  on  algebraic  steps  rather  than  ones  involving 


* 
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differential  equations.  For  algebraic  operations  there  Is  sane  hope  of  computer 
mechanization  so  that  the  generality  criteria  of  section  1 could  then  be  met; 
but  where  each  step  still  Involves  partial  differential  equations  there  Is  none. 

(c)  The  obvious  assumption  consistent  with  (a)  and  (b)  above  Is 

<p  = f (A)  + f^A)  cos  u + f2(x)  COS  2w  + ...  (A20) 

since  <t>  must  clearly  be  an  even  periodic  function.  However,  substituting 
this  Into  (A13)  and  separating  the  harmonics  leads  to  an  Infinite  sequence  of 
equations  each  of  which  contains  linear  differential  operators  of  second  order. 
Solution  of  such  a system  in  practical  terms  seems  out  of  the  question,  and 
Its  Intractable  nature  arises  from  the  strong  coupling  between  equations  — in 
the  sense  that  the  expression  for  f contains  all  f's  for  other  values  of 
n.  This  In  turn  Is  a consequence  of  the  addition  formulas  for  trigonometric 
products  and  so  cannot  be  avoided  by  any  assumption  of  the  form  (A20). 

(d)  It  appears  therefore  that  (A20)  is  precisely  the  wrong  kind  of 
assumption,  and  of  course,  the  coupling  problem  arising  from  (A20)  can  be 
alleviated  by  using  Chebyshev*  polynomials  instead  of  trigonometric  series, 
since  then  the  major  algebraic  steps  reduce  to  manipulations  with  power  series. 
Putting 

C = cos  (j,  n - eA 

then,  the  analog  of  (A20)  Is 

<P  = 9Q(n)  + g-j(n)^  + g2(n)?2  + ...  (A21) 

which  greatly  reduces  the  coupling  problem  but  doesn't  solve  It.  In  fact,  use 
of  (A21)  leads  to  a system  In  which  the  equation  for  gn  involves  gn+2 
which  Is  again  undesirable  since  each  differential  equation  cannot  be  solved 
outright  In  general  terms  which  might  conceivably  lead  to  a reduction  to 
algebraic  steps. 


*In  a quick  spot-survey  of  16  mathematics  textbooks  and  assorted  dictionaries, 
the  author  found  8 different  spellings.  This  was  the  most  popular  being  quoted 
on  37.5%  of  occasions! 
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Examining  the  reduction  to  a suitable  s,n  form  shows  what  the  proper 
procedure  should  be.  Noting  that  y has  factor  sin  w,  (A12),  let  Y be 
defined  as 


where 


and  u„(0 


y = sin  <„  • Y 


B B 

Y ( 5 » n ) = n + ~ + — ( ? ) + ... 
n 

is  the  Chebyshev  polynomial  of  second  kind,  viz, 

U (r)  = + 1 

n'*'  sin  a) 


(A22) 

(A23) 


Using  (A22)  and  (A23)  in  (A13)  then  gives  the  transformed  differential  equation 

(A24) 


Y (V2,  - & ) + n2Y  * + (1  - c?)Y 

n n s c. 


0 


with  the  boundary  conditions 


= 0 when  n = 1 


<j>  ^ £n  when  D -►  or> 


(A25) 


But  if  (A21 ) is  applied  to  (A24)  it  can  immediately  be  seen  that  the  offending 

factor  which  produces  a shifting  upward  of  powers  is  the  term  (1  — i 
. „ 2 SS 

coming  from  transformation  of  the  v <*>.  On  the  other  hand  the  n-part, 

n <J>nn  + n<J>fi  preserves  the  same  powers  in  each  differential  operator. 


A5.  THE  ASSUMPTION  FOR  * 


(a)  On  the  basis  of  section  A4d  above,  the  assumption 

4>  = nF0(0  + ^ F^c)  + \ F2(C)  + ...  (A26) 

n 

is  more  logical  than  (A21). although  both  (A20)  and  (A21)  appeal  more  immedi- 
ately to  intuition.  Using  (A26)  it  is  possible  to  produce  a system  which  can 
be  reduced  entirely  to  programmable  algebraic  steps.  This  reduction  is 
elementary  but  involved,  and  only  the  important  stages  can  be  recorded  here. 


Putting  (A26)  progressively  into  (A24)  leads  to  the  first  step  that 

V2$  - = n40  + ^ A1  + ...  -J-  An  + ... 

n 


where  An  Is  the  differential  expression 

4n  5 (1  - E2)F»  - 2Ef;  + n2F„  , n > 1 

The  temptation  to  use  this  as  a basis  for  Iteration  (as  opposed  to  the  use  of 
the  contribution  from  v <p  alone)  on  the  grounds  that  It  may  introduce  non- 
two-dimensional  effects  quicker  If  one  divided  (A24)  by  Y,  can  be  resisted 
very  easily  by  noting  that  the  first  member  would  be 

40  = (1  - C2)F"  - 2r,F'  + F0  » 0 (A27) 

which  does  not  have  a solution  proportional  to  e,  as  it  must  if  (A25)  is  to 
be  satisfied. 


Continuing  the  development,  we  note  instead  that  the  product  term  involv- 
ing v24>  can  be  reduced  in  the  following  way.  Write 


[ ®2^1  j C Pi  3p  j 

Y = njl  + -j  + — j—  + ...(  = nil  + — + -Tj-  + . . . I - n(l  + p) 

* n n \ ‘ n nJ  1 

and 


2 ( A1  A2  ) 

7 ♦ ~ ■ njAo  + T + T-  **•  [ = n(A0  + A) 

say,  so  that  the  product  becomes 


Y {v2<>  - e*  } = n2(l  + S)(a0  + a)  = n2(l  + 8A0  + A + SA) 

where  the  final  term  pa  is  0(n"^).  The  other  two  terms  in  (A24)  yield 
similar  expressions  with  final  products  of  0(n~4);  in  detail,  if 


y 


n 


♦n 


= 1 


2s2 

7 “1“  ** 


n n 


F,  2F2 

~!~~r 


then 


i2Yn«On  = n2(l  - &HF0  - h = n2  {Fq  - ^Fq  - F + Sf  } 
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and,  if 


then 


ej 

T+  ~3  + •• 

n n 

Fj  FJ 

Fi  + 4 + T + 

n n 


“ nB' 


= n(F;  + F') 


(1  - C2)Y^c  - n2(l  - S2)(B’F'  + 3 1 F * } 

2 

so  that  ultimately  a factor  n can  be  cancelled  and  the  differential  equa- 
tion reduced  to 


{Aq  + 0AO  + A + Fo  -$Q  -?  + (1  - ?2)0'F'}  + [f3A  + BF  + (1  — 52)b'F']  - 0 

(A28) 


which  can  be  written  for  convenience  of  discussion  as 


{H}  + [G]  = 0 (A29) 

In  this  expression  the  collection  G contains  only  terms  of  at  least  &(n”4) 
and  the  collection  H starts  from  order  zero.  This  separation  of  the  higher- 
order  products  turns  out  to  be  a very  important  step,  since  isolated  study  of 
the  dependence  of  G on  the  functions  Fn  shows  how  the  whole  solution  can 
be  reduced  to  mechanizable  algebraic  steps. 


(b)  However,  the  structure  of  G can  be  appreciated  only  by  first  examin- 
ing the  behavior  of  H.  Interpreting  H from  (A27)  and  the  definitions  of 
a,b,F»F  shows  that 

H H 

h = (a0  + F0)  + 4 + "f  + •••  (A3o) 

n 

where 

Hn  * 4n  - "F„  + Vn  ~ Fo"Bn  + <’  ~ <«’ > 

Since  G is  0(n~^)  it  follows  that  the  solutions  for  the  first  three  Fn's 
are  independent  of  G and  are  given  by 

HQ  =0,  H1  = 0,  H2  = 0 (A32) 
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and  in  particular  the  first  of  these  according  to  (A29)  and  the  definition  of 
aq  (A27)  yields 

0 - S2)F JJ  - 2^  + 2F0  = 0 

This  is  a Legendre  equation  of  degree  1 with  the  solutions 

F0  = AP] (c)  + BQ1 (c) 

where  P-j  ( 5 ) and  (e ) are  the  Legrendre  polynomials  of  degree  1 and  A,  B 

ore  arbitrary  constants.  In  fact, 

P-|(c)  = 5»  Q-|(c)=^-Cln  j — 1 

and  since  the  second  kind  [Q-|(c)]  gives  a logarithmic  singularity  at  £ = +1 
we  set  B - 0 on  the  grounds  that  4>  must  be  a smooth  function.  Then 
choosing  A = 1 gives 

Fo  = A0  = (*33) 

so  that  the  first  term  of  (A26)  becomes  n5  which  generates  a proper  free 
stream  as  n °°»  as  it  should  according  to  (A25). 

Substituting  this  back  Into  (A31 ) shows  that  Hn  can  be  reduced  to 

Hn  ' An  - nFn  ~ (n  + ’>«„  + (1  " 

” Ln  * ,),;Un-l(5)  “ (1 

where  Ln  is  the  Legendre  differential  operator  of  degree  (n  - 1), 

Ln  - (1  - «2)F"  -2«f;  + n(n-  l)Fn  (A34) 

Furthermore,  the  expression  { } Is  a well-known  identity 
(n  + lkU^U)  - (1  - C2)U^_1  (^)  = nUn(f.) 

so  that,  finally 
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and  therefore  the  complete  sequence  of  solutions  from  (A29)  must  satisfy 

L1  - B1U1 


L2  = 2B2U2 
L3  " 3B3U3  " G3 


(A35) 


where  we  have  written 


L„  = nBr.Un  ~ G„ 

n n n n 


G3  G4 

G = T + T + 

n n 


since,  as  already  noted,  6 is  0(ti  ). 


(c)  With  regard  to  the  Gn  quantities  arising  from  the  higher-order 
products  it  can  be  seen  from  (A28)  that 

G = BA  + 8F  + (1  — £Z)B'F' 

and  each  of  these  products  can  be  expanded  into  a power  series  in  1/n  by 
using  the  definitions  of  b,a  etc.  to  give 


W =Z>r  <VlVl-r  + r<"  + ’ - ^r-lW  + <’  - «VlFJ+l-rJ  <A36> 

r=l 


The  details  of  how  this  can  be  simplified  sufficiently  to  enable  algecralc 
computation  to  be  affected  for  arbitrary  n will  be  considered  later  (section  9). 
For  the  moment  the  important  feature  is  that  Gn  depends  only  on  Fn2,  An_2, 
etc.,  so  the  only  coupling  is  to  functions  already  known  from  previous  lower 
orders  — which  is,  of  course,  the  vital  point  of  contrast  to  all  the  other 
Iterative/analytical  methods  discussed  in  Section  4. 


Now  the  general  solution  of  Ln  = 0 which  remains  finite  at  c = + 1 is 


“nPn-l<5> 


S3 


where  an  Is  an  arbitrary  constant.  Furthermore,  a particular  solution  of 

Ln  ’ "W‘> 

Is 

F„  “ 

since 

(1  -52)TJ-25T;  + n(n-l)Tn=  -nUn 
where  Tn(s)  Is  the  Chebyshev  polynomial  of  first  kind 

Tn(e)  = cos  nu 

Thus  F^F2,/y\2  are  polynomials  In  e and  so,  therefore.  Is  G3  according  to 
(A36),  However,  It  Is  not  necessarily  obvious  that  every  F„  Is  thereby  a 
polynomial  since  F3  Involves  the  particular  solution  for  a right  side  (G3) 
which  Is  (essentially)  an  arbitrary  polynomial  rather  than  a special  firm. 

But  the  highest  degree  In  G3  is  1 so  we  could  represent  the  problem  for 
F3  as 

L 3 e 3 B 3 U 3 p|  • 3 Af:  + B (say) 

to  emphasize  this  aspect.  A particular  solution  for 

t-3  = P-t . F3  = q1  (say) 

is  q-j  = ( 1 /4 ) Ac  + ( 1/6 )B  which  Is  still  of  degree  1 only.  The  highest  degree 
In  G4  Is  c2,  but  now  we  know  that  the  highest  degree  In  Gg  is  still 
only  because  the  solution  for  F3  did  not  Increase  the  degree  of  p-j. 

Proceeding  In  this  manner  It  Is  clear  that  _1f  a particular  solution  qn_2 
can  be  found  for  the  equation 

Ln  * Pn-2 

where  pn_2  Is  at  most  of  degree  £ then  every  FR  will  be  a polynomial 
of  degree  <n. 
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(d)  To  see  that  this  Is  the  case,  consider  the  superposition  of  solutions 
of 

Lf|(q)  » ^ (A36) 

where  q Is  assumed  to  be  a power  series 

q - q1  + q2s  + q2<-2  + ...  (A37) 

which  may  or  may  not  terminate*.  Operating  on  (A37)  with  Ln>  and  equating 
5 coefficients  gives  the  system  of  equations 

Nq-|  + 2.1q3  « 0 | 

(N  - 2.1 )q2  + 3.2q4  = 0 I 

> N - n(n  - 1)  (A 38) 

(N  - m(m  - 1 ) + m(m  + 1 )qm+2  = * I 

• o • I 

whose  solution  depends  on  the  value  of  m In  relation  to  n. 

If  m » n then  q^  can  be  found  from  (A38)  and  therefore  all  the 
higher  coefficients  displaced  by  two.  There  Is  no  way  to  find  the  coefficients 
displaced  by  one,  nor  those  for  m,  m-2,  m-4,  ....  However,  since  we  are  only 
Interested  In  finding  an£  special  solution  these  can  be  taken  as  zero  and  the 
recursive  computation  of  q^,  q^g  ...  provides  a satisfactory  procedure  — 
except  for  the  fact  that  this  sequence  does  not  terminate. 

Evidently  If  m > n all  coefficients  up  to  m + 2 must  again  be  taken 
as  zero  with  the  same  consequences. 

Fortunately,  the  case  for  m < n is  different  since  q^2  = 0 and  we 
may  take  all  powers  n,  n-2,  ...  as  zero  until  m Is  reached.  Therefore, 
the  recursive  sequence  starts  with  q^  and  descends,  viz. 


*The  choice  of  numbering  In  the  coefficients  may  appear  a little  unusual,  but 
this  particular  system  will  be  used  frequently  since  It  Is  consistent  with 
the  computer  numbering  and  numerical  steps. 
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until  or  q2  is  reached. 

Thus  it  is  possible  to  obtain  a special  solution  to  (A36)  which  is  a poly- 
nomial of  degree  m - 1 provided  i«i  does  not  exceed  n - 1.  As  can  be  seen 
from  section  cf  this  is  just  the  case  required  where  the  polynomial  on  the  right 
(pn2)  is  of  degree  at  most  n - 2. 


Therefore,  we  have  established  Inductively  that: 

1.  Gn  = pn_2  is  of  degree  at  most  n — 2 

2.  qn_,,  is  of  degree  at  most  n - 2 

3.  Fn  is  always  a polynomial  of  degree  at  most  n 

On  this  basis  it  is  more  convenient  to  use  the  p,q  notation  and  write  the 
system  of  equations  (A35)  as 


L 

L 

L 


1 

2 

3 


B,U, 

2B2U2 


3B3U3  " P1 


(A40) 


The  solution  of  this  system  is, 
discussion. 


L - nB  U - p o 
n n n rn-2 

• • • 

according  to  section  c,  and  the  above 


f,  * »,P„<0 
F 2 = U2P1  ^ ^ BgT  g ( C ) 

F3  = 2^^  B3T  3 ^ — 


f = n,r  (f)  — B T > r 1 — n 


(A41 ) 


All  of  these  are  elementary  polynomials  and  the  only  items  standing  in  the  way 
of  a general  automatic  computation  scheme  are  the  details  of  the  pn,  qp 
calculations  and  the  manner  In  which  the  boundary  conditions  must  be  satisfied 
In  order  to  determine  the  constants  «n. 

Both  of  these  items  together  with  the  general  behavior  of  the  solution 
can  be  clarified  by  a brief  examination  of  the  first  few  approximations. 

A6.  STRUCTURE  OF  THE  FIRST  FEW  TERMS 


(a)  Having  seen  that  4>Q  = nFQ  gives  a free  stream,  it  is  instruc- 
tive to  see  next  what  the  inclusion  of  the  next  two  terms  yields  since  these 
are  still  independent  of  any  p -*  q operations.  Let  «-|  denote  the  approxi 
mation  up  to  F^  so  that 

4>i  = nC  + (a-j  — BiO  = n COS  <d  + ~ (a-j  — B.j  cos  u) 

from  (A41 ) . Then 


-a.|  + (1  + B-j ) cos  u 


and  it  is  clear  that  there  is  no  choice  of  o-j  which  makes  this  vanish  for 
all  h>  unless  it  happens  that  B-|  = -1,  Bn  = 0(n  >•  1).  in  that  case  if 
a-j  = 0 the  mapping  represents  a zero  thickness  (plate)  profile  parallel  to 
the  stream  direction. 


Taking  the  next  term 

^2  = + o’  ^al  _ B-jTj(s)]  + [agC  “ 8212(c)] 


• - a,  ♦ 
n l 


/ B1  a2  \ B2 
^ n — + — >5- j cos  hi  — — cos 


2h) 


so  that 


-»|  + (1  + B1  - Zo^)  cos  ui  — 2B2  cos  2u> 
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It  Is  Impossible  to  choose  and  a2  to  make  this  vanish  for  all 
unless  B2  * 0.  If  Is  zero,  then  the  choice 


1 + B, 


gives 


/ B1  1+Bl\ 

V ” + 17”/C0S  w 


The  case  of  ellipsoids  Is  considered  In  more  detail  In  the  next  section,  but 
at  this  point  It  Is  worth  noting  that  If  every  Bn  = 0 (Including  B^ ) the 
mapping  generates  a circular  profile  and  under  these  circumstances 


(nti7)c°s“‘(r  + i7)' 


which  Is  the  correct  potential  for  a sphere. 


For  any  other  body  $2  Is  an  approximation.  It  is  not  correct  for  an 
ellipsoid  when  Bj  f 0 even  though  every  other  Bn  Is  zero  for  such  a body. 
It  Is  this  curious  feature  that  makes  the  ellipsoids  worth  further  study  In 
their  own  right.  If  B2  t 0 then  $2  Is  a poor  approximation  In  general 
because  the  boundary  condition  will  be  In  error  In  the  cos  2W  term. 

(b)  Since  it  is  not  difficult  to  work  out  p^  and  p2  by  hand  we  next 
consider  Leaving  out  the  details,  the  calculated  values  of  p-j  and  p2 
obtained  from  (A36)  and  , F2  are: 

P]  = ^B-ja'j,  P2  = (48^2  + 6B2aj )? 

and  then  the  procedure  of  section  5d  gives 

Blal  1 

d]  = “j ' » ^2  = F ^2Bla2  + 3B2a1  ^ 

so  that  the  whole  solution  for  <j>4  Is 


= + 


i («,P0  - B,T,)  * \ («,2P,  - B2T2)  ♦ \ (o3P2  - B3T3  - 

n n 

\ (a4P3  “ B4T4  " F C2Bla2  + 3B2u1^) 

n ' 


Putting  the  Legendre  polynomials  In  terms  of  cos  nu  and  applying  the  boundary 
condition  then  gives  for  each  frequency  In  decreasing  order 


1 


— | a4  + 3B3  - ° 

— j + 2 EL,  * 0 

— ^ (28^2  + 3B2<^-| ) “ 2a2  + Bj  + 1 B 0 
“T  a3  + °1B1  “ “1  = 0 


(A44) 


with,  as  before,  the  understanding  that  the  cos  4u>  term  cannot  be  accounted 
for.  The  first  two  equations  give  Immediate  solutions 


!B3’ 


8 R 

tb2 


and  then  the  last  two  reduce  to  the  2x2  matrix 


(8B1  - 10)a2  + 128^  = 9B3  - 5(1  + B] ) 


3(B|  - l)Ql  » 2B2 


which  can  be  solved  as 


2 B2 


“I  T 


1 8B2 
“2  = (5B'”-”TUT  9B3  * 


- 5(1  ♦ B 


"1 


(A45) 


This  shows  that  the  <*n  coefficients  will  not  In  general  depend  linearly  on 
the  Bn  coefficients, and  that  some  kind  of  small  matrix  solution  will  be 
required  for  higher-order  terms.  In  fact,  the  algebra  up  to  the  4x4 
matrix  solution  was  carried  out  by  hand,  but  only  a few  key  points  are  recorded 
for  use  In  the  next  section. 


The  structure  of  <j»g  Itself  In  terms  of  Fn,  Is  clear  enough  from  the 
basic  formula  (A26)  together  with  (A41).  Expressions  for  q-j  and  q2  have 
already  been  given  and  the  results  for  q3  and  q4  follow  from  the  calcula- 
tion of  p3  and  p4  through  (A36)  and  the  procedure  of  section  5d.  The 
results  are: 
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1 


1 


where 


13  ■ Co  + C2«*«  <<4  ■ °1E  * D3'3 


Co  * “ TO  [30Bla3  - 4B2«2  + (MB1  + <4B3>01  ] 
C2  * yjy  + + 1 ^ 


(A46) 


°1  -"W 


56Bi(i4  + 30B? 


>a3  + ("T^  B1  + 24B3)a2 


(’ 


4 (100B4  4 


irBiB?)“i 

Oj  3 yjy  [20B-ja4  + + 24B3a2  + 40B4o-j] 

Satisfaction  of  the  boundary  condition  on  the  surface  leads  to  six  equations 
for  the  six  o*n  constants  which  again  break  down  into  two  explicit  formulas 
for  ac  and  a*,  viz. 


_ 5 x 128  n 

*6  ~ rrsr  Bs» 


4 x 64 


a5  = 5 x 35  B4 
and  4 simultaneous  equations  for  o^,  a2»  “3*  “4 

«3(45B1  - 63)  + a2(50B2)  + ^(SOB-j)  3 64B4  - 56B2 

a3(75B1  - 105)  + a2(270B2)  + ^ (140  - 140B1  - 180B3  + 140B2)  3 -144B4 


(447 ) 


a4(60B1  -90)  + a3(63B2)  + a2(72B3>  * “iC120^)  K lOOBg  - 108B3 


(448) 


a4(420B-j  — 63)  + a3(1305B2)  + a2^43  + a1^44  = l^OOBg  ~ 420(1  + B-j) 
where  the  two  coefficients  C43  and  C44  are 

C43  3 (-840  4 672B1  4 1800B3  - 648B2) 

C44  3 (1008B2  4 1210B4  - 1872B-jB2) 

These  clumsy  results  are  Included  only  for  the  algebraic  arguments  of  the  next 
section  and  for  completness  In  the  sense  that  they  are  invaluable  guides  in 
working  out  and  checking  general  automatic  computing  schemes. 
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A7.  COMPARISON  WITH  ELLIPSOID  SOLUTIONS 


(a)  The  reason  that  It  Is  worthwhile  considering  ellipsoids  In  detail  Is. 

of  course,  that,  they  represent  the  only  simple  general  class  of  parametrically 
varying  axl symmetric  bodies  with  exact  solutions.  Furthermore,  there  Is  a 
curious  feature  Involved  when  comparing  the  classical  exact  solutions  with  the 
sequence  •••  developed  In  the  foregoing  sections.  This  arises  In  the 

following  way.  With  all  BR  * 0 except  Bj»  separation  of  real  and  Imagin- 
ary parts  In  (A9)  and  putting  x ■ 0 gives 

x$  ■ (1  — Bj)  cos  w,  y$  - (1  + B1)  sin  u 

If  we  choose  to  take  C * 0 and  so  shift  the  origin  to  the  usual  central 
point  for  ellipses.  Thus,  keeping  the  x-axIs  always  aligned  with  the  flow, 

semi -major  axis  = a = 1 - 

semi-minor  axis  * b a 1 + B-j 

showing  that  -1  £ Bj  <. +1,  with  B-j  » 0 giving  a circular  profile  as  noted 
In  section  6b,  and  *1  i Bj  <_  0 giving  the  usual  range  of  thickness  ratios 
representing  "ovary"  ellipsoids.  Now  It  Is  apparent  that  the  foregoing  theory 
gives  •••  directly  In  terms  of  B-j  Irrespective  of  whether  B1  Is 

> or  <0.  However,  the  classical  results  (e.g.,  Mllne-Thomson  [A7])  show 
that  the  analytical  form  of  solution  Is  different  according  as  a > b or 
a < b,  l.e. , according  as  B-j  < 0 or  Bj  > 0. 

This  situation  can  be  clarified  by  "expanding"  the  classical  solutions  In 
some  way  and  showing  that  the  expansions  so  obtained  agree  with  ‘P04>2^4  ••• 
Irrespective  of  the  sign  of  B^. 

(b)  The  classical  solutions  can  be  expressed  by  means  of  elliptic  coord- 
inates as  follows: 

o When  a > b (-1  1 Bj  <_  0)  put  B^  =*  -e"2lJ  and  write  y * cosh  (x  + y), 
yo  = cosh  u then 

♦ ■ 2«'u “ jr - i/j j 

(A49) 
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r 


o 


When  a < b (1  >_  >_  0)  put  Bj  ■ e'*^  and  write  o - slnh  (x  + p), 

aQ  ■ slnh  u then 


4>  ■ Ze“v  cos  w 


la  tan~^(l/o)  - 1] 

[tan*  (1/oJ  ” ort/(on  + !)] 


(A50) 


'o'  'o'  ' 0 

A reasonable  supposition  Is  that  If  (A49)  and  (A50)  are  capable  of 


yielding  expansions  which  can  be  matched  to  ♦ , $2»  $4  •••»  then  the  case  of 
large  x should  be  the  one  to  examine.  The  first  term  of  (A49)  Is 


e~u  (e*+n  * cos  u * (ex  + e“2pe"x)  cos  w 


cos  u 


In  the  earlier  notation.  Thus,  temporarily  denoting  the  constant  In  (A49)  by 
h we  have 


cos 


If  x and  p are  both  taken  as  large  then  x + u Is  large  and  It  Is  natural 
to  put  e = e~^x+V1^  and  find  the  expansion  for  small  e.  The  term  In  x + p 
Is  then  expressible  as 


Y y + 1 , _ 1 /l  . \ /I  ♦ e\  ,.42.  8 4 . 12  6 . 

*7\7+  e)  1n\r^r)"l  * Ie  *T5E  +55Te  + — 

.2  = --2p«-2a 


or,  since  e = e 


-Bj/n2. 


cos 


f / Bl\  4.  h /4Bl  1 8B^  1 + lZBl  1 \l 

“Ir  " r) + h - TT- 7 + ir  7 - -)J 


(A51) 


In  order  to  derive  consistent  comparisons  for  Individual  orders  from  (A51)  It 
Is  convenient  to  proceed  In  a slightly  unusual  manner  because  of  the  appar- 
ently very  complicated  nonlinear  dependence  of  the  4>^g  terms  on  B1 
through  the  <*n  constants  (see  e.g.  (A43)  and  (A44)  for  and  o^). 

Since  h is  the  constant  appropriate  to  an  Infinite  number  of  terms.  It  seems 
reasonable  to  argue  that  one  should  designate  h^*  h^,  hg  as  the  constants 
appropriate  to  the  orders  n-2,  n~*»  n”6  and  then  these  can  be  determined  by 
noting  that  every  order  must  Individually  satisfy  the  surface  boundary 
condition.  This  procedure  Is  certainly  consistent  with  the  philosophy  under- 
lying the  structure  of 


i 


( 
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For  Instance,  on  this  basis  (A51)  yields 


$2  a cos 


“ [("  - r) 


4B.  . 
n 


(A52) 


and  then  we  must  have 


3(1  + Bj ) 

— w, — 


(A53) 


Similarly  ^ and  $g  are  given  by  (A51)  with 

-15(1  + B.| ) 
h4 = ; 


h6  = 


105(1  + B.) 

y-  (A54) 

8B-|  (35  - 28B]  + 27B*) 


For  the  other  case  (A50)  let  the  constant  be  denoted  by  k then 


4>  = cos  a) 


and  again  let  e = e"^A+p^  so  that 


(i-m- 


k(o  tan”1  - — 1 ) 


(° tan''  7_,)  = 7^~£) Un''  r^r-1 

i — e 


'b  — 


4 

3 


8 

IS 


E — 


12 

3* 


This  time 


2 

E 


2A 


and  If  this  Is  substituted  Into  the  expansion  above  the  same  expression  as 
(A51 ) results.  Since  the  constants  k2»  k^,  kg  ...  would  be  chosen  by  the 
same  argument  It  follows  that  (A51)  through  (A54)  represent  the  successive 
orders  of  approximation  for  both  cases,  thus  resolving  one  aspect  of  the 
curious  feature  mentioned  In  section  a. 


(c)  Next  It  Is  necessary  to  show  that  the  <f>2,  $g  of  section  6 are 
consistent  with  (A51)  through  (A54).  In  the  case  of  ^ It  Is  obvious  by 
Inspection  on  comparing  (A42)  with  (A52).  For  put  Bg.  B3  to  zero  and 
then  «3  and  are  zero  according  to  (A44).  Furthermore,  according  to 
(A45)  aj  = 0 and 

5(1  + Bj ) 

“2  = “ - 10) 
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so  that  If  these  results  are  substituted  into  (A43)  the  formula  for 
becomes 


*4 


(A55) 


which  can  be  seen  to  be  identical  to  (A51)  with  h^  given  by  (A54). 


Finally,  for  <p g with  all  Bn  coefficients  zero  except  B-j,  the  equa- 
tions (A47)  and  (A48)  for  the  an's  reduce  to 

cig  = 0,  ag  = 0,  a^=Q,  a2=0,  a-|  - 0 

and 

105(1  + B-| ) 

2 (210  - 1688^  + 162B2) 

Then  the  expression  for  <j>g  becomes  quite  simple  because  q3  and  q^ 
simplify  to 

q3  = 0,  ^4  = ~ If  B2a2c 

thus  giving 

a6  = cos  w [(n  " 7r)  + “2  ( T “ F B1  V + 3T  B1  V)]  (A57) 

Using  the  «2  given  by  (A56)  in  (A57)  it  is  easy  to  see  that  this  result 
agrees  with  (A51)  and  (A54)  showing  that  the  classical  expansions  agree  with 
the  algebraically  derived  approximations  <t>2»  of  this  theory  for 

ellipsoidal  bodies. 


(d)  Finally,  the  comparisons  will  be  complete  if  we  can  show  that  the 
constants  h2>  h^  and  hg  are  progressi vely  better  approximations  to  h 
and  k,  the  exact  constants  appropriate  to  an  infinite  number  of  terms. 
Denoting  e”u  by  c (small),  the  first  form  of  constant  from  (A49)  is 
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h = 


2e 


-u 


H/2  In  [(y0  + D/(y0  - 1)]  - C(y0)/(yo  ~ 1)]> 


which  can  be  written  In  the  forms 

1 


h = 


1 

^[8/3  + (24/5) e2  + (48/7) e4  ...] 


0 + Bj ) 


(A58) 


B-j [8/3  - (24/5)B1  + (48/7)Bj  - ...]  Bj[8/3  - (32/15)B-,  + (72/35)Bj  ...] 

(A59) 

2 2 

since  c = - -B-j.  Progressively  including  higher-power  terms  in  [ ] 

gives 


3(1  + B-j ) 


’ 


15(1  + B-j ) 105(1  + Bj) 

» ~ 5 

'1'  D 8Bj(35  - 28Bj  + 27Bj) 

2 . 


which  are  identical  to  (A53)  and  (A54).  The  expansion  for  k in  e is 

the  same  as  (A58)  except  for  alternating  signs  which  give  exactly  (A59)  again 
? 

since  this  time  e = +Bj. 


Thus  both  forms  of  classical  formula  give  the  same  expansion  for  the 
constants  of  each  type.  Of  course  the  two  forms  for  the  classical  formulas 
result  from  the  same  general  argument  in  more  general  complex  variables  so 
none  of  this  Is  really  very  surprising.  Nevertheless  this  topic  cannot  be 
discussed  further  here.  We  have  shown  that  ellipsoidal  bodies  are  correctly 
represented  by  the  sequence  of  potentials  in  this  theory.  It  is  more 
appropriate  to  inquire  about  the  convergence  of  this  sequence,  and  this  can 
best  be  discussed  in  terms  of  the  velocity. 

A8.  FORMS  OF  THE  VELOCITY 

(a)  The  general  theory  for  $ as  given  in  sectionsAl  through  A7  will  be 
of  little  value  unless  we  can  be  sure  that  the  velocity  — and  in  particular 
the  surface  velocity  — can  be  computed  accurately  — without  problems  from  the 
trailing  edge  singularity  for  instance.  For  the  present  it  is  only  the 
surface  velocity  (a  = 0)  which  is  of  Interest,  and  denoting  magnitude  by  Q 
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we  have  from  section  3 


q((U)  ■ - h^tbt  (f^.c 

In  the  case  of  the  conformal  mapping  used  here 

jdzl 

fa! 


h = h .1^-1 
nl  "2 


io) 


and  on  the  surface  c = e so  that 

dz  f 

h!<0)  ■ 3 ETI 


_ ds 
<Jw 


where  s Is  the  arc  length  along  the  profile  measured  from  the  trailing  edge 
(which  is  taken  as  the  origin  in  all  general  computations).  Thus  using  s as 
a subscript  for  surface  value 

3<f>  /3a> 

% - ~ -&7G  (A60) 


and  the  need  for  caution  arises  from  possible  zero  values  of  ds/dw  as  a 
consequence  of  the  assuned  forms  for  the  mapping  derivative  [(A7)  and  (A8)]. 

Setting  c = e1u)  In  (A7)  and  (A8)  gives 

t = tt  (blunt)  = £(1  + ip  cos  + . ..)2  + (&2  sin  2u>  + . ,.)2j 

(A61) 

= 0 (cusp)  = 2 sin  j £(1  + cos  u>  + a2  cos  2u>  + ...)2 

2ll/2 

+ (sin  u + a2  sin  2u  + J (A62) 

so  that  ds/dtu  does  Indeed  vanish  at  w = 0 when  t = 0.  Ji  Is  convenient 
to  note  that  (A61 ) and  (A62)  can  be  expressed  In  terms  of  the  surface  modulus 
of  the  g(c)  function  of  section  3 for  both  cases,  that  Is: 


Sb = M“>l  1f  T s " 

= 2 sin  2 |9s(w)|  if  t = 0 


(A63) 


66 


where  I g$ (oi)  | Is  a positive  periodic  function  which  has  no  zeroes  or 
Infinities  - as  Is  obvious  from  (A61)  and  (A6 1). 

Now  the  form  of  <j>  given  by  (A 26)  and  (A41)  has  the  natural  terms  con- 
sisting of  cos  nos  and  powers  of  cs  '.os  <o.  Thes.  latter  can  always  be 
expressed  In  terms  of  cos  w,  cos  2<o  ...  by  using  the  reverse  transformation 
for  the  Chebyshev  polynomial  of  first  kind,  so  that  It  Is  always  possible  to 
express  <f>  In  the  form 

<f)  a ( n ) + G'j(n)  cos  to  + Ggfn)  cos  2a>  + ... 

by  a simple  reorganization.  Consequently, 

|^  = — ^(n)  sin  a)  + 2G^ ( n ) sin  2o>  + ...J 

from  which  a factor  sin  w ^an  be  entracted  leaving  a modulation  term  which 
Is  a function  of  the  Chebyshev  polynomials  of  second  kind.  That  is, 

* — s1nu»G(c,n)  (A64) 

o u) 

where 

G(4,n)  = G-|  ( n ) + 262(0)11,(5)  + 3G3(n)U2(c)  + ...  (A65) 

so  that  the  velocity  from  (A60),  (A63),  (A64)  can  be  written 

Qs  ■ 1f  T ■ ” 

• " - 0 (c.«P) 

s 

and  both  of  these  forms  are  computationally  trouble  free  for  all  u>.  As 
expected  blunt  bodies  have  stagnation  points  at  both  ends,  whereas  a body  with 
a cusp  trailing  edge  has  a finite  nonzero  velocity  at  that  trailing  edge 

(ui  * 0). 
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(b)  For  the  ellipsoids  the  results  of  section  7 can  be  used  to  give  some 
Idea  of  the  maximum  velocity  convergence.  It  Is  Interesting  to  compare  the 
successive  terms  with  the  two-dimensional  value  which  can  be  denoted  as  the 
zero-th  order  approximation  very  conveniently.  The  complete  solution  for  two 
dimensions  Is, of  course, 

4>0  = (n  + cos  w = 2 cosh  a cos  <*> 

so  that,  dropping  the  now  unnecessary  "s“  subscript,  the  zero-th  order 
approximation  to  surface  velocity  is 


n _ 2 sin  ai 

Qo  ■ TFJ^TT 


or 


(A65) 


Irrespective  of  the  shape.  For  the  order  2,  4,  6 terms  on  ellipsoids  the 
G function  surface  values  follow  from  (A42),  (A55),  and  (A57)  leading  to 
the  sequence 


G 


o 


2 


(1  — B-g  ) 
(1  — B1) 

(1  ~ B] ) 


+ i (1  + B,) 

*5 


, / 35  - 14B,  + 9B \ \ 

1 0 + B])  ( j) 

c 1 \35  - 28B1  + 278^  / 


(A66) 


As  an  illustration,  the  case  of  ellipsoids  of  fineness  ratios  a/b  = 3/1  and 
a/b  = 1/3  are  given  in  the  table  below.  Exact  values  for  G are  computed 
from  (A491  and  (A50)  by  setting  A = 0 to  give  the  slightly  simplified  forms: 


G = 


4e" 


2 cosh  u - sinhu  In  [(cosh  v + l)/(cosh  u - 1)] 


i < 0 - -e"£‘u(a  > b) 


= n t — ; B,  > 0 - e"2p  (a  < b) 

cosh  y tan  (1/sinh  p)  — sinh  p 

whilst  the  values  of  | gs (oj)  | at  w * t r/2  corresponding  to  the  point  of 
maximum  velocity  are: 


li 
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when  a = j , b = |-  then  B-j  * + tj-  ; Jg$(j)|  =»  ^ 

The  results  are  summarized  In  the  table  below  with  the  % errors  defined  as 

* (approx  value  - exact  value)  x 100 

% error  - exacTl^Ui L 


Table*  of  G Values  and  Maximum  Velocities  for  Various  Orders 
of  Approximation  on  Two  Ellipsoids 


Slender  Ellipsoid  a/b  = 3 

Blunt  Ellipsoid  a/b 

= 1/3 

Order 

G values 

vfe 

%°Error 

G Values 

vfc 

% §BSr 

0 

2.0 

1.333333 

+18.8 

2.0 

4.0 

+45.8 

2 

1.75 

1.166667 

+4.0 

1.25 

2.5 

-8.8 

4 

1.714286 

1.142860 

+1.9 

1.50 

3.0 

+9.4 

6 

1.698431 

1.132300 

+0.9 

1.317568 

2.635136 

-3.9 

Exact 

1.682953 
, ..  ■ 

1.121970 

1.371324 

2.742648 

As  an  illustration  not  too  much  should  be  deduced  from  these  tabular 
values.  Perhaps  we  should  just  remark  here  that  the  signs  oscillate  for  the 
excessively  blunt  case  (a/b  = 1/3),  but  not  for  the  more  usual  type  of  body 
(a/b  a 3/1),  and  that  convergence  is  adequate  but  not  Impressive.  In 
section  9 detailed  results  will  be  presented  for  a variety  of  bodies  so 
judgment  should  be  deferred  until  then. 

(c)  Before  passing  onto  general  computation  one  other  remark  concerning 
accuracy  Is  worth  drawing  attention  to.  Naturally  the  convergence  of  velocity 
is  crucially  dependent  on  the  convergence  of  the  Bn  coefficients.  Some  idea 
of  the  behavior  of  these  mmbers  can  be  gained  by  brief  consideration  of  the  < 

Joukowski  profile  case  on  which  Kaplan  (section  4)  obtained  such  good  answers 
with  a scheme  essentially  sensitive  only  to  83.  A Joukowski  profile  has  the  * 

exact  modulating  function 
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($ee[A4])  where  d is  real  for  a symmetric  case  and  0 < d < 1.  Therefore, 
the  expansion  coefficients  are 

an  - ndn_1  - (n  - l)dn;  Bn  - -d"’^!  - d)2 

When  | d | small  It  Is  clear  that  the  convergence  of  the  Bp  numbers  is 
good.  For  Instance  in  the  case  quoted  by  Kaplan  d * 0.15,  giving 

B]  - -0.7225,  B2  * -0.1084,  B3  ■ -0.01626,  B4  = -0.00244 

The  dominant  "ellipsoid-effect"  or  leading  term  Is  proportional  to 
1 + Bj  = 0.2775  so  the  amplitude  ratio  of  the  first  neglected  term  in  Kaplan’s 
calculation  would  be  essentially  proportional  to 

0.00244  1 

0775“  < m 

which  goes  some  way  toward  explaining  the  satisfactory  results  obtained  with 
such  limited  sensitivity  to  the  Infinite  sequence  of  Bn  coefficients. 

However,  It  Is  not  difficult  to  get  cases  of  practical  Importance  in 
which  the  Bn  convergence  Is  very  much  worse  than  this, as  will  be  seen  in 
section  10,  where  test  cases  for  the  general  solution  having  deliberately 
chosen  peculiarities  are  presented. 

A9.  GENERAL  CALCULATION  SCHEME 

(a)  The  overall  structure  of  the  general  calculation  scheme  Is  based  on 
the  sequence  of  operations  already  described  algebraically  In  the  earlier 
sections  (5  through  6).  The  computer  program  was  written  In  Fortran  IV  for 
use  on  the  Douglas  IBM  370/168  system  and  can  be  regarded  as  proceeding  by 
the  following  stages: 

1.  Input  Bn  coefficients  and  other  data  necessary. 

2.  Compute  double  subscripted  arrays  for  the  coefficients  of  special 
polynomial  used  (Pn(c),  Tn(e), 
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3.  Fill  the  first  three  orders  F^,  Fg,  Fj  of  the  solution  variables 
using  the  algebraic  answers  discussed  In  Section  5b  and  c.  All 
solution  polynomials,  unlike  the  special  polynomials  In  2,  are 
functions  of  the  ap  constants  and  so  must  be  represented  by 
3 subscript  arrays  which  store  the  coefficients.  Thus,  F,  p,  q,  ... 
are  all  typically  represented  as 


Fn<F->  - F„.l«> 


+ F„>2(E)a,  + 


* Fn,rvH(5>“n 


where  each  polynomial  Fn  ^ ( C ) Is  of  the  form 

F j ( 5 ) = f_  j i ^ f_,  ^ oC  ^ ...  f_  j »ii £ 

n,1  n,1,l  n,l,z  n»i,n+i 


so  that  the  coefficient  array  f ^ k requires  storage  N,  N + 1, 

N + 1 If  N Is  the  highest  order  desired.  Since  three  subscripted 
arrays  use  a lot  of  core  space, some  trouble  was  taken  In  the  program 
to  minimize  the  number  of  such  arrays  actually  used.  The  current 
version  only  requires  two;  but  even  so  a great  deal  more  could  be 
done  to  improve  the  efficiency  of  this  particular  version  of 
program. 


4.  Starting  with  n = 4 each  Fp  coefficient  array  Is  computed  using 
(A41)  and  the  p and  q functions  of  order  up  to  n — 2.  In 
addition  the  p function  of  next  highest  order  is  computed  at 
this  stage  and  the  coefficients  stored  for  future  use. 

5.  A subroutine  computes  q coefficients  from  p coefficients  by 
using  the  algorithm  of  section  5d  and  these  are  also  stored  for 
future  use  as  In  4. 


6.  When  the  Fn  functions  have  been  computed  up  to  the  order  desired 
(N  say),  a separate  subroutine  handles  satisfaction  of  the  boundary 
conditions  by  first  converting  fn  ^ k to  ?n  ^ k such  that 

F„  *(?)  * i t + f 2 oT,  + f„  i ,T,  + ...  f * _.,T 

which  Is  an  appreciable  simplification  at  this  stage  since  Tk  is 
merely  cos  k«.  Using  the  arguments  of  section  6,  a set  of  N — 2 
simultaneous  equatlnns  for  the  ap  constants  follows,  which  can  be 
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solved  by  a simple  standard  matrix  Inverse  of  low  order  operating 
wholly  within  core. 

7.  A final  subroutine  computes  and  outputs  the  surface  potential  and 
velocities  using  the  formulas  of  section  8. 

This  description  of  the  essential  steps  followed  by  the  program  is  only 
Intended  as  an  outline  for  discussion  purposes;  but  most  of  the  programming 
involved  in  the  ommitted  details  is  elementary  and  does  not  deserve  comment. 

However,  one  aspect  of  this  scheme  does  deserve  brief  amplification. 

(b)  The  central  operation  which  includes  all  the  "difficult"  part  of  the 
whole  theory  (essentially  all  the  nonlinear  aspects)  is  the  computation  of  p 
from  the  extant  Fn  functions  according  to  (A36).  It  has  already  been 
remarked  that  not  only  does  the  whole  theory  hinge  on  the  fact  that  each  F 
depends  (via  p)  only  on  previously  computed  F's,  but  that  also  it  is 
possible  to  prove  that  FnU)  is  a polynomial  of  degree  at  most  n 
(section  5c,  d).  However,  the  boundary  condition  procedure  can  always  be 
reduced  to  a matrix  inverse  of  order  (N  - 2)  x (N  — 2)  with  and 
determined  directly,  because,  in  fact,  the  pp_2  which  effects  Fn  is 
actually  only  of  order  n - 3 rather  than  n - 2. 

In  section  5c  it  was  noted  that  if  the  special  solutions  q were  poly- 
nomials given  a polynomial  p of  degree  at  most  n - 2,  then  the  validity 
of  the  whole  procedure  could  be  proved  Inductively.  That  Pn(-Sn+2^  itself 
is  of  degree  at  most  n is  obvious  from  (A36),  but  the  additional  vanishing 
of  the  highest  order  term  requires  a little  closer  examination.  Denoting 
n — 2 by  m,  (A36)  can  be  rewritten 
m 

pm  Br  {Ur-lVl-r  + + 1 ” r)Ur-lFm+l-r  + (1  ” 2jUr-l Fm+l-r} 

r=1  (A67) 

and  each  major  contribution  is  a product  of  two  polynomials.  A special  poly- 
nomial product  subroutine  keeps  track  of  the  coefficients  for  this  frequently 
used  operation,  but  for  the  highest  order  terms  alone  the  answer  can  be 
written  down  directly.  Nevertheless  (A67)  is  not  the  best  form  of  either  the 
numerical  construction  or  the  algebraic  scrutiny.  Using  the  expression  for 
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(section  5a),  It  Is  easy  to  reorganize  (A67)  Into 


A_  In  terms  of  F„ 
n n 

m 

pm  ■ Z)  sr  [a?  0 - * <m  * »(■  ♦ 1 - r)V,  Wr] 

r=l 


(A68) 


which  is  the  form  actually  used  In  the  program.  Furthermore,  it  is  obvious 
from  the  structure  of  Fn  (A41 ) that  the  highest  order  term  is  always 
contributed  by  -BnTn(?)  and  so  is  always 

Consequently  the  second  part  of  (A68)  contributes  to  c™  the  amount 


m 

BrBm+l-r(m  + + 1 - r)2m~]  (A69) 

r=l 

The  highest  order  term  In  U iF^j  is,  by  the  same  reasoning, 

-Wr<m  * 1 - r)2m-1en-1 

so  that  the  contribution  to  cm  from  the  first  part  of  (A68)  must  be 


m 

+ BrBnH-l-r^m  + + 1 “ 

r=l 

which  exactly  cancels  (A69),  showing  that  the  sm  term  in  pm  vanishes. 
Thus  the  computing  scheme  can,  as  explained,  make  use  of  this  Information  to 
keep  the  an~matrix  as  small  as  possible.  The  cases  of  N = 4 and  6 as 
worked  out  in  section  6b  are  illustrations  of  this  feature. 


A10.  SOME  EXAMPLES  OF  DIRECT  VELOCITY  CALCULATIONS 

(a)  In  order  to  test  the  theory  a nunber  of  profiles  were  generated  by 
using  a modification  of  the  "pole  airfoil"  theory  reported  in[A4]»  This  pro- 
cedure is  extremely  flexible  since  poles  and  zeros  of  the  mapping  derivative 
can  be  chosen  at  will  within  the  unit  circle  in  order  to  generate  bumps  and 
dimples  on  the  profile.  The  point  of  all  this  is,  of  course,  to  provide  test 
cases  of  sufficient  character  to  reveal  more  information  about  the  accuracy 
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than  can  be  obtained,  for  Instance,  by  only  using  Joukowskl  profiles.  We  have 
already  seen  that  this  can  be  very  misleading  (section  8c)  and  since  most 
common  test  shapes  are  members  of  this  larger  class  of  pole  airfoils  (profiles) 
It  Is  considered  to  be  a useful  adjunct  to  the  whole  theory. 


For  the  test  cases  described  below  the  form  of  dz/dc  started  out  as  a 
quotient  of  isolated  linear  groups,  viz. 


dZ_  (t-C ,)(t-l2)  ...  U -cB) 

3t  - U--d,)U-J2J  ...  (t-dj 


|c„l.  |<l„l  < 1 (A70) 


where  the  first  factor  can  be  used  to  generate  blunt  or  cusped  shapes  - since 
the  latter  is  the  only  allowable  example  of  the  singularity  being  on  the 
surface;  the  choice  then  being 


c-|  = 1 , dj  = 0 

to  give  the  usual  type  of  cusp. 


Otherwise  various  criteria  apply  to  the  choice  of  Cp  c2»  . ... 
d^,  d2,  ...  required  to  generate  specific  effects  which  need  not  be  explained 
here  in  detail  (see[A4]),  except  to  note  that  the  generation  of  extreme  surface 
convolutions  can  only  be  achieved  In  general  by  having  some  |cn|  and  some 
| dnJ  very  close  to  unity. 


Equation  (A70)  can  be  expressed  as 


ds  . 
Hi" 


dki  < i 


(A71 ) 


which  is  the  fundamental  pole  form  where  the  quantities  x(  are  the  pole 
strengths.  A simple  partial  fraction  program  determined  these  strengths  from 
(A70),  and  then  (A71)  can  be  used  to  give  an  exact  formula  for  the  shape.  In 
addition  (A71)  lends  Itself  to  an  Immediate  expansion  which  Is  uniformly 
convergent  in  |c|  >.  1 since  every  |d|<|  < 1,  viz. 


dz  _ 
Hi" 


(A72) 


< 


< 
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where 


m 


k-1 


(A73) 


These  coefficients  show  that  the  ultimate  rate  of  convergence  Is  controlled  by 
the  largest  |dk|  and  therefore  It  follows  from  the  remarks  above  that 
extreme  convolutions  will  have  poor  convergence.  In  the  context  of  this  work, 
of  course,  everything  depends  on  the  significance  of  the  words  "extreme"  and 
"poor."  This  subject  Is  best  left  to  the  examples  given  below,  except  to  note 
in  passing  that  the  Bp  coefficients  used  for  the  theory  here  are  defined 
directly  from  (A 72)  by  integration.  Thus 


z 


C + 


agreeing  with  (A9),  and  so 


(remembering  that  a^  * 0 for  closure)  which  shows  that  the  Bn  coefficients 
converge  faster  than  the  an  coefficients. 


It  should  be  noted  that  this  system  of  test  cases  generates  exact  shapes 
and  exact  BR  coefficients  whose  rate  of  convergence  can  be  controlled  and 
studied.  Clearly  this  Is  necessary  for  the  testing  phase,  but  It  does  not 
represent  a bias  In  favor  of  the  current  method.  Any  profile  enclosing  a 
schllct  singly-connected  domain  can  be  mapped  Into  the  unit  clrcld.and  It  Is 
a trivial  matter  to  write  such  a mapping  program  and  to  make  the  associated 
Bn  coefficients  available  for  input  to  this  theory.  In  this  sense  neither 
is  the  special  system  used  here  a restriction  of  generality.  A mapping  pro- 
gram is  available  for  finding  these  coefficients,  but  was  not  used  in  the 
cases  described  below.  This  question  will  be  considered  further  in  the  dis- 
cussion of  the  Inverse  (design)  mode  (section  11). 


b.l.  The  first  figure  (Fig.Al)  shows  what  might  be  described  as  a very 
bland  basic  case  (body  1)  in  which  the  profile  was  a smooth  blunt 
shape  generated  by 
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c. |  » 0.9, 

d,  = 0, 


c2  * -0.5 
d2  - 0.4 


The  velocity  distribution  based  on  only  10  Bn  coefficients  is  shown 
compared  with  the  very  accurate  "higher-order  Neumann"  method  of 
Hess  [A2 ) . Evidently,  agreement  is  good  for  a trivial  Investment  of 
computer  effort  because  the  "hardest"  task  involved  in  this  theory 
for  only  10  terms  is  the  inversion  of  an  8x8  matrix*.  For  the 
Neumann  run^  and  output  stations  of  the  method  of  this  appendix,  101 
points  on  the  half  body  were  used.  Of  course,  the  blandness  of  this 
case  leads  one  to  suspect  that  it  Is  particularly  favorable  to  the 
kind  of  theory  advocated  here,  just  as  the  Joukowski  case  used  by 
Kaplan  (see  section  8c)  turned  out  to  be  particularly  favorable  in 
that  context.  In  fact,  the  rate  of  convergence  for  this  case  ^s 
unduly  good.  Atthe64th  station  (very  close  to  the  maximum  velocity) 
the  zero-th  order  approximation  (two-dimensional)  is  1.520655  and  the 
succeeding  orders  are: 


Order 

2 

4 

6 

8 

10 


9s 

1.226915 

1.227366 

1.225472 

1.225760 

1.225739 


showing  nrecisely  the  kind  of  weighted  choice  one  should  avoid  in 
test  cases 1 


2.  However,  the  results  obtained  for  the  second  body  (Fig.  2)  continue 
to  show  excellent  agreement  even  though  only  10  terms  are  used.  This 
body  is  a perturbation  of  the  basic  case  (1)  In  which  a bump  was  gen- 
erated deliberately  near  the  nose  - as  can  be  seen  from  Fig.  A2. 

The  resultant  sharp  peak  In  velocity  is  produced  accurately  in  detail 
and  appears  to  be  a good  demonstration  of  the  excellent  convergence 
properties  of  the  method. 


♦There  is  little  point  in  trying  to  isolate  proper  computer  CPU  times  from  these 
runs  because  they  are  so  small  that  the  figures  are  largely  obscured  by  vari- 
ous accounting  and  other  semloptional  features.  For  Instance,  In  this  run  the 
CPU  time  was  quoted  as  0.011  minutes  and  the  1/0  as  0.044,  but  velocities  of 
order  0,  2,  4,  6,  8 and  10  were  individually  computed  and  stored  which  is 
typical  of  many  options  not  always  needed. 
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3,  Body  No.  3 Is  a further  perturbation  of  even  more  extreme  character 
(Fig.  A3).  As  can  be  seen,  the  velocity  shows  quite  violent  response, 
to  the  artificially  created  bump-dip-bump,  but  the  theory  again  gives 
remarkably  good  results  with  only  10  terms  when  compared  with  the 
higher-order  Neumann.  Obviously,  this  case  is  a still  better  demon- 
stration of  the  relative  insensitivity  of  the  basic  method  to  higher 
order  Bn  convergence. 

4.  The  final  case  is  interesting  because  it  shows  how  a cusp  is  handled 
by  the  theory.  The  body  is  essentially  similar  to  case  2 except  that 
the  Bn  sequence  terminates  at  Bg  and  that  there  is  a cusp.  It 
can  be  seen  (Fig.  A4)  that  agreement  (again  with  10  terms)  is  excel- 
lent everywhere  except  close  to  the  trailing  edge.  A much  greater 
number  of  terms  Is  required  to  give  a good  answer  in  the  Immediate 
locality  of  the  cusp.  This  is  exactly  what  would  be  expected  on  the 
grounds  of  the  discussion  given  at  the  end  of  section  2 where  it  was 
pointed  out  that  if  the  structure  of  the  Bn  coefficients  gave  good 
answers  for  the  nonsingular  (blunt)  case  then  obviously  it  would  be 
less  satisfactory  for  cases  where  a further  zero  term  should  be 
extracted. 

These  cases  cannot  be  considered  exhaustive  and  evidently  further  work 
should  be  done,  particularly  for  arbitrary  bodies  where  the  mapping  program 
to  determine  the  Bn  coefficients  Is  required  first.  To  some  extent  the 
design  cases  considered  next  relieve  any  residual  anxiety  that  perhaps  the 
above  results  were  subtly  (unintentionally)  biased  In  favor  of  this  theory. 
Accepted  at  face  value,  however,  the  above  results  show  very  clearly  the 
excellent  convergence  and  economy  of  the  theory. 

All.  DESIGN  (INVERSE)  OPERATION 

(a)  As  pointed  out  in  Section  1,  one  of  the  important  advantages  of 
analytical  solutions  In  engineering  is  the  readiness  with  which  they  can 
usually  be  converted  to  a design  (inverse)  mode.  In  the  context  of  airfoil 
theory  it  happens  that  the  design  problem,  where  a desired  velocity  is  given 
and  an  associated  shape  determined,  is  both  exact  mathematically  and  non- 
iterative. Airfoils  designed  in  this  way  are  those  which  - in  a least 
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squares  sense  - give  velocities  "closest"  to  the  desired  Input  and  as  such 
have  recently  proven  to  be  of  great  practical  value  [A8]. 

In  fact,  interest  in  design  methods  for  axisymmetric  flow  was  one  of  the 
original  stimuli  for  this  study,  particularly  as  there  was  no  design  procedure 
of  any  kind  available  at  the  start  of  this  work  except  In  two  dimensions. 
However,  recently  a procedure  for  designing  axisymmetric  shapes  by  simply 
iterating  a direct  (Neumann  type)  method  has  finally  been  perfected  by  Bristow 
[A9]  and  we  are  fortunate  to  be  able  to  compare  a few  cases  of  Bristow's  method 
and  the  procedure  advocated  here. 

But  first  a brief  description  of  the  Inverse  mode  for  the  current  theory 
is  needed.  Once  again  theimportant  feature  is  the  decoupling  of  the  purely 
geometric  aspect  and  the  velocity  aspect.  In  two  dimensions  this  decoupling  is 
complete,  whereas  for  ax i symmetric  flow  the  gradient  of  <f>  still  depends  on 
the  Bn  coefficients.  Thus,  unlike  the  two-dimensional  case,  the  design 
procedure  for  axisymmetric  flow  is  Iterative.  The  steps  used  in  the  current 
program  version  are  very  similar  in  character  to  those  needed  to  determine 
the  Bn  coefficients  (i.e.,  the  mapping)  from  a given  profile,  or  to  find  the 
flow  about  a given  (two-dimensional)  airfoil.  As  such  they  constitute  an 
iterative  scheme  belonging  to  a general  class  whose  convergence  and  accuracy 
properties  have  excellent  precedents  from  two-dimensional  experience.  Briefly 
the  procedure  is  based  upon  the  central  Idea  that  if  ds/du>  can  be  found  at 
a certain  stage  then  the  appropriate  an  coefficients  can  be  calculated 
directly  by  employing  a transformation  which  is  common  in  two-dimensional 
mapping  theory.  Looking  at  (A61)  and  (A62)  it  might  be  thought  that,  given 
ds/dw,  there  is  no  way  to  determine  the  an  coefficients  directly,  but  this 
is  not  the  case.  The  argument  by  which  these  coefficients  can  be  determined 
is  typical  of  the  great  power  and  flexibility  of  mapping  methods  and  is  rarely 
appreciated  as  such. 

(b)  Briefly,  if  dz/de  Is  given  by  series  of  the  forms  (A7)  and  (A8), 

then 

x = tt  (blunt)  In  ||  = In  gU)  = g(c) 

t » 0 (cusp)  In  * In  (!—■£■)  + In  g(e)  * In  (1  - 1)  + f(e) 
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and  in  both  cases  g (?)  Is  an  analytic  function  for  |c|  >_  1 which  -*•  1 as 
|c|  ->  ».  Thus  f(c)  is  an  analytic  function  for  |?|  ^1  which  +0  as 
|cj  ■*  00  and  can  be  represented  by 


In 


g(0  = In  (l  + ^-+  ...)  = f(c)  = ~ + ^|  + 


(A74) 


In  a relation  of  this  form  the  determination  of  the  ap  coefficients  given 
the  bn  or  vice  versa  is  a purely  algebraic  procedure  which  is  easy  to 
compute.  On  the  other  hand,  the  connection  between  the  surface  values  of  f 
and  g is  f = In  g so  that  if  f = u + iv  (say)  then 

o j d j j 


“s  = ln  9S  ! 


»s  * ar9  9S 


Therefore,  using  (A63),  the  first  of  these  gives 


u = ln  Si 
s doj 


(t  = tt) 


■ ar- ■»"  (**<»?)  i'-v 

can  be  calculated.  But  according  to  (A74)  u_  is  the  real 

^ 5 

part  of  the  bn  series  with  real  coefficients  (for  the  ax i symmetric  case  — 
by  symmetry)  giving 


from  which  u 


us  = b-j  cos  (d  + b2  cos  2u  + ...  (A75) 

where  b-j  = a^  * 0 if  t = tt  or  b-j  = a-j  = 1 if  x = 0 due  to  the  closure 
condition  (A5). 

Thus  it  follows  that  the  bn  coefficients  can  be  determined  directly 
by  Fourier  analysis  of  u$  and  therefore  the  ap  coefficients  indirectly  by 
inversion  of  (A74)  as  remarked  before.  In  addition,  the  magnitude  of  b-j  is 
clearly  a measure  of  "compatibility"  at  any  star  of  the  iteration  or,  in  the 
particular  case  of  two-dimensional  design,  a direct  measure  of  the  plausibility 
of  the  desired  velocity  distribution.* 


♦The  inverse  problem  is  not  unique  in  either  two-dimensional  or  axisymmetric 
flow.  The  best  one  can  expect  is  some  minimun  measure  of  closeness  to  the 
desired  input  and  a least  squares  Fourier  criterion  has  been  Found  very  satis- 
factory in  two  dimensions  (see,  e.g.  [A8]. 


(c)  The  steps  in  the  Iterative  design  mode  for  the  theory  developed  here 
can  then  be  described  In  the  following  terms: 


Step  1.  Guess  a set  of  Bp  coefficients.  The  runs  described  below  all 


use  those  for  body  1 (section  lOd)  as  a starting  system. 


Step  2.  Calculate  3$s/9 to  by  using  (A64)  and  (A65).  For  this  step  the 


major  part  of  the  basic  direct  program  as  described  in  section  9a 
is  used  as  a subroutine.  Only  number  7 is  left  off. 


Step  3.  Use  (A60)  to  give  ds/du  and  hence  usU). 


Step  4.  Find  the  bn  coefficients  from  (A75)  by  Fourier  analysis  of 


us(w)  according  to  section  b above.  For  this  step  a small  Fast 


Fourier  transform  subroutine  tailored  for  ordinary  engineering  usage 
in  Fortran  IV  was  used.  (In  a modern  computer  context  this  has  to  be 
recognized  as  one  of  the  big  advantages  of  Fourier  methods  in 
general ). 


Step  5.  From  the  bn  coefficients  calculate  the  ap  coefficients  and 


hence  the  Bn  coefficients  which  in  turn  define  the  new  velocity 


and  associated  profile. 


Limits  on  convergence  were  defined  by  testing  the  RMS  changes  in  y^.,  and 


Q..  Experience  has  shown  that  usually  velocity  is  the  more  sensitive  param- 

J 


eter  and  that  it  is  very  important  not  to  impose  a too  stringent  requirement 
on  accuracy.  In  fact,  a limit  of  0.001  is  quite  sufficient  for  most  engineer- 
ing purposes  and  is  usually  safe.  The  importance  of  adequate  caution  ibr  itera- 
tive procedures  of  this  kind  derives  from  the  fact  that  such  schemes  very  often 
contain  numerical  steps  of  a very  commonplace  nature  (e.g.,  interpolations,  quad- 
ratures, ...)  which  are  in  fact  much  less  accurate  than  the  overall  theory  or 
truncation  thereof;  but  which  frequently  get  overlooked.  This  remark  does  not, 
cf  course,  apply  to  the  Fast  Fourier  Transform  routine  used  here  which  returns 
11  significant  digits.  However,  this  particular  theory  does  have  a nonlinear 
dependence  on  the  Bp  coefficients,  no  scaling  on  the  higher-ordar  contributions 
(where  the  coefficients  of  the  Legendre  and  Chebyshev  polynomials  become  large) 
and  a very  early  truncation.  These  are  reasons  enough  for  not  overdoing  the 
convergence  criterion. 
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A12.  RESULTS  OF  INVERSE  CALCULATIONS 


(a)  As  an  obvious  first  check  the  Inverse  mode  was  run  on  all  the  direct 
cases  described  before  (section  A10).  Convergence  was  remarkably  good  again 
with  only  10  terms,  and  most  of  the  results  could  not  be  distinguished  from 
the  Input  on  any  reasonably  plottable  scales.  Body  3 was  the  only  one  with 
sufficient  character  to  be  interesting  and  it  was  the  first  one  to  be  run 
with  direct  comparison  with  the  method  of  Bristow  [A9].  Hence  it  is  the  only 
one  worth  more  discussion  here. 

In  fact.  Fig.  5 shows  the  original  exact  body  and  exact  velocity  (higher- 
order  Neumann)  used  for  input  to  the  design  process.  It  also  shows  the  results 
of  13  cycles  of  the  current  thoery  compared  with  20  cycles  of  Bristow's  method. 
From  this  figure  it  appears  that  after  20  cycles  the  Bristow  method  still  has 
a long  way  to  go  and  that  quite  small  deviations  In  velocity  near  the  nose  are 
associated  with  quite  large  errors  In  overall  thickness. 

On  the  other  hand,  the  method  advocated  here  Is  very  good  after  only  13 
cycles  (with  no  under-relaxation  factor,  see  section  b below).  Furthermore, 
it  is  quicker  since  the  CPU  times  were  approximately  in  the  ratio  1:10. 

(b)  Another  case  of  design  operation  is  shown  in  Fig.A6.  It  is  particularly 
interesting  in  that  It  represents  a more  typical  situation  in  which  the  designer 
inputs  the  very  crude  straight  line  distribution  (Fig.  A6)  knowing  full  well 
that  it  is  impossible  in  detail,  but  that  It  represents  the  kind  of  levels  and 
gradients  he  requires.  This  figure  shows  that  after  7 cycles  of  the  current 
theory  and  20  cycles  of  Bristow's  the  desired  Input  is  very  nearly  being 
achieved.  However,  it  is  clear  that  the  two  body  shapes  do  differ  near  the 
trailing  edge.  This  is  a reflection  of  the  same  situation  as  encountered  in 
section  a above  - namely  that  the  Bristow  method  approaches  more  slowly  to 

the  limit,  a fact  that  is  made  apparent  on  Fig.  A6  by  noting  the  closeness  of 
the  agreement  between  the  higher-order  Neuamnn  answers  for  the  7 cycle  body  and 
the  7-cycle  velocities.  Once  again  the  CPU  times  were  approximately  in  the 
ratio  1:10. 
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Some  details  about  this  run  are  worth  noting.  Firstly,  it  should  not  be 
thought  that  such  a close  match  to  the  desired  Input  is  possible  In  general. 

The  case  shown  Is  one  which  was  known  to  be  close  to  a "real"  velocity  distri- 
bution by  experimentation.  Had  the  desired  Initial  rise  and  trailing  edge 
values  been  associated  with  a level  of  say  1.2  Instead  of  1.05,  then  the  design 
theory  would  have  shown  how  far  from  possible  such  a request  was.  It  would  also 
have  returned  some  "closest"  match,  leaving  it  up  to  the  designer  to  use  this 
information  to  Introduce  a compromise.  Secondly,  the  remarkably  good  converg- 
ence shown  for  the  current  theory  may  not  always  be  achieved  in  practice. 

Fig.  A7  shows  the  behavior  of  the  B-j  coefficient  and  RMS  y error  as  functions 
of  cycle  nunber  with  and  without  an  under-relaxation  parameter  of  1/2  whose 
need  was  rather  obviously  Indicated.  There  are  some  Intuitive  grounds  for 
regarding  an  "average"  as  a good  choice;  but,  of  course,  there  is  no  mathemat- 
ical ground  for  assuming  that  this  value  is  optimum  In  any  rigorous  sense. 
Obviously  more  experience  with  the  method  is  desirable.  Thirdly,  it  should  be 
noted  that  by  ordinary  engineering  (fluid  mechanic)  standards,  these  straight 
line  Inputs  are  excessively  crude  having  discontinuities  In  slope  and  other 
features  quite  uncharacteristic  of  analytic  functions.  It  Is  therefore  a very 
convincing  argument  In  favor  of  mapping  and  Fourier  analysis  methods  that  they 
can  treat  such  unfavorable  distributions  smoothly  and  accurately. 


A13.  CONCLUSIONS 


V 


4 


A theory  for  the  general  solution  of  axisymmetric  Incompressible  potential 
flow  in  terms  of  analytic  functions  of  elementary  type  has  been  presented.  For 
the  cases  usually  of  Importance  to  engineering  (blunt  and  cusp  trailing  edges) 
comparison  with  existing  numerical  methods  has  shown  that  with  only  10  terms 
excellent  results  can  be  obtained  for  negligible  expenditure  of  computer  time 
on  bodies  with  more  character  than  usually  encountered.  When  used  in  the 
design  (inverse)  mode  the  method  is  apparently  quite  satisfactory  being  both 
quicker  and  more  accurate  than  the  iterative-Neuamnn  algorithm  of  Bristow 
which  is  (as  far  as  Is  known)  the  only  alternative  nonlinear  design  process 
for  axisymmetric  flow. 

Clearly  more  work  needs  to  be  done  to  explore  the  full  range  of  this  method, 
and,  maybe  It  would  be  worthwhile  studying  the  solutions  for  arbitrary 
(0  < t < -..)  trailing  edge  angles.  The  results  so  far  are  certainly  encouraging. 
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extreme  deformation  of  the  basic  shape  showing  how  the  complicated  velocity  distribution 
n well  matched  with  only  10  terms  when  compared  with  the  higher  order  Neunann  method. 


A body  generated  by  a terminating  B-sequence  with  a cusp  trailing  edge  showing  the  analytical 
solution  (10  terms)  compared  with  the  higher  order  Neumann  method. 


Bristow  (20  cycles). 


APPENDIX  B 


A NUMERICAL  STUDY  OF  SOME  OPTIMUM  DRAG  CONFIGURATIONS  IN  TWO  DIMENSIONS 

B1 . INTRODUCTION 

As  a contribution  to  the  overall  investigation  of  minimum  drag  bodies  of 
revolution,  certain  two-dimensional  studies  were  undertaken.  The  main  purpose 
of  these  was  to  gain  Insight  into  the  possibilities  Inherent  in  the  mapping 
approach  to  profile  definition  which,  it  was  felt,  might  yield  a sufficient 
injection  of  analyticity  into  the  whole  problem  for  a genuine  variational 
analysis  to  become  possible.  As  a technique  for  separating  geometrical  and 
flow  aspects,  conformal  mapping  is  of  great  importance  in  two  dimensions,  and 
its  value  in  axisymmetric  flow  as  a simplifying  mechanism  can  be  judged  by 
Appendix  A where  an  analytic  method  for  flow  analysis  or  design  cf  bodies  of 
revolution  is  developed  ab  initio  by  using  unit  circle  mapping.  At  very  least 
this  kind  of  mapping  constrains  the  allowable  form  of  velocity  in  such  a way  as 
to  focus  attention  on  relatively  few  parameters  characterizing  "mode  functions" 
of  the  correct  type. 


Although  it  would  seem  that  this  is  a necessary  prerequisite  for  vitroductlon 
of  calculus  of  variations,  it  was  found  that  in  the  time  available,  no  concrete 
variational  approach  could  be  formulated  for  the  axisymnetric  case.  However, 
the  two-dimensional  pilot  studies  did,  indeed,  lead  to  a number  of  interesting 
results  and  advances  In  understanding  which  can  be  summarized  as  follows: 

a.  A two-dimensional  complete  variational  formulation  does  appear  to  be 
possible.  However,  development  of  the  underlying  mathematical  tech- 
niques is  insufficient  for  such  a problem  to  be  solved  at  the  present 
time.  Variational  methods  using  complex  variable  function  theory 
appear  to  be  both  much  more  subtle  than  the  conventional  real  variable 
methods  and  totally  unexplored.  We  shall  not,  therefore,  discuss  this 
subject  further  here. 

b.  It  is  possible  to  reduce  both  the  axisymmetric  and  two-dimensional 
problems  to  minimization  problems  involving  a definite  sequence  of 
parameters  to  be  chosen  such  as  to  minimize  certain  combinations  of 
integrals. 


PI 


c.  A very  simple  computer  program  was  written  to  find  solutions  of  such 
minimization  problems,  and  the  two-dimensional  applications  of  this 
program  under  various  restricted  circumstances  led  to  some  interesting 
results  - not  the  least  of  which  was  the  inherent  feasibility  of  this 
approach  for  the  unique  definition  of  the  minimum  drag  shape  given 
enclosed  area  or  volume. 

The  work  described  in  (b)  and  (c)  was  not  completed  owing  to  insufficient  time, 
but  the  feasibility  is  adequately  illustrated  by  the  results  described  below. 

B2.  TWO  DIMENSIONAL  DRAG  INTEGRAL 

According  to  Thwaites[Bl]  as  interpreted  by  A.M.O.  Smith,  the  drag  of  a 
two-dimensional  symmetric  shape  in  wholly  turbulent  flow  can  be  expressed  by 
Spence's  integral  as 


Drag  _ 2 

(1/2 ) pu£c 


0.0242SR"1/5 


(Bl) 


where  p is  the  half  perimeter,  c the  chord  and  Q the  surface  velocity 
divided  by  expressed  as  a function  of  arc  length  s.  In  this  formula  the 
Reynolds  number  R is  based  on  chord  so 

R = Uwc/v  (B2 ) 

where  Uro,  v have  their  usual  significance.  (Bl)  can  aTso  be  written 


2(0.0452R“1/6 


Q4ds 


5/6 


(B3 ) 


where  the  expression  in  ( ) is  the  Blasius  formula  for  a flat  plate  and  the  non 
dimensional  group  in  [ ] is  a convenient  parameter  worth  a name,  viz. 

c 

'o  * If  <M> 

0 
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In  making  comparative  geometric  studies  it  is  clear  that  both  the  chord  and 
free-stream  velocity  should  be  kept  constant,  but(B3)  alone  is  not  very  useful 
as  it  stands,  since  Iq  >_  1 so  that  minimization  of  Cqc  will  always  yield 
the  flat  plate  solution.  Alternatives  which  could  conceivably  be  useful  from  a 
practical  engineering  standpoint  would  keep  (say)  enclosed  area,  maximum  thick- 
ness or  perimeter  constant  as  well.  Some  experience  with  these  has  shown  that 
the  area  constraint  produces  the  most  desirable  type  of  optima,  and  the  remainder 
of  this  study  is  concerned  only  with  drag  optima  for  constant  area,  chord  and 
free  stream. 

As  far  as  the  conversion  of  this  to  a readily  computed  form  in  terms  of  the 
mapping  theory  is  concerned,  the  area  factor  must  be  nondimens ional i zed  in  order 
that  body  size  should  not  be  an  influence,  so  a sensible  parameter  for  minimiza- 
tion is 


C°A  = CDc 


c 

SK 


where  A = total  enclosed  area 


But  the  quantity  most  easily  computed  is  the  half  area  integral  which  can  be 
readily  expressed  by 


and  then 


c 

“T  " "T  f ydx  = *a 

2cZ  cZ  J A 


(say) 


giving  from  (B3)  and  (B4) 


Cqa  = (0.0452R'1'/6* 
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Finally  for  reasons  of  partial  comya^fsons  w^Ch  otH«r  r>;ioUs  it  y»s,  decided  that 
the  integrals  based  on  the  C9nfigyr;w*f*  J>e  used,  which  would  giva 
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CDa  = (0.0452R‘1/6)  21/6(2Id)5/6/(2Ia)7 


with  the  purely  geometric  expression 


F = 


(21ft) 
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separable  as  a measure  of  configuration  efficiency.  In  terms  of  F then 

CDa  = (0.0452R"1/6)  21/6F  (B8) 


and  the  work  described  below  is  concerned  with  the  minimization  of  F. 


B3.  MAPPING  FORMAT 

In  Appendix  A a fairly  complete  description  of  the  mapping  format  required 
for  transformation  of  the  profile  in  the  z-plane  to  the  unit  circle  in  the 
5-plane  is  given.  Although  the  basic  philosophy  concerning  constraints  on  this 
mapping  function  is  unchanged  for  the  application  here,  an  appreciable  generali- 
zation is  desirable.  Thus  the  work  of  Appendix  A was  restricted  to  circumstances 
under  which  it  was  feasible  (possible)  to  work  out  complete  solutions  for  general 
axlsymmetrlc  flow,  and  in  fact  only  blunt  noses  with  either  blunt  or  cusped 
trailing  edges  were  actually  analyzed. 

As  far  as  the  numerical  analysis  of  the  integrals  of  section  B2are  concerned, 
however,  it  is  sensible  to  begin  by  statinq  the  most  general  case  — that  for  which 
finite  angles  and  Tg  are  permitted  at  both  ends  of  the  body.  Under  these 
circumstances  the  mapping  derivative  can  be  written 

£-(1 t7)U2[’4t7+-]  <»> 

where 

P-j  ~ 1 — Tj/.,  Mg  = 1 — 12>/ir 


and 
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In  order  to  ensure  closure.  In  addition  It  can  be  shown  by  symmetry  that  for 
profiles  representing  a body  of,  revolution  every  an  coefficient  must  be  real 
and,  of  course,  the  surface  behavior  of  the  derivative  Is  obtained  by  setting 
c - e^  for  points  on  the  unit,  circle.  In  fact  the  connection  between  the 
mapping  derivative  on  the  surface  and  the  arc  length  is  of  central  importance, 
and  follows  from  the  observation  that.  If  subscript  "s"  denotes  surface  value, 
then 


(BIO) 


(Bll) 


where  s is  the  profile  arc  length  measured  from  the  trailing  edge,  and  eg 
is  the  geometric  surface  angle  (or  flow  angle). 


There  are  a number  of  different  formulas  which  can  be  used  for  the  various 
integrals  arising  in  section 82  which  will  be  very  useful  in  deriving  certain 
special  and  exact  cases.  However,  there  is  very  little  that  can  be  done  with 
the  fundamental  integral 


c 


unlike  the  chord  and  area.  Thus  for  this  integral  we  note  that  the  surface 
velocity  (in  two  dimensions  only)  can  be  written 


n 2 sin  w 
Q ' 117^ 


3U) 


so  that 


I, 


TT 

-If 


2 sin 
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D c ^ (ds/dw) 
Putting  c = e1u)  in  (B9)  and  (Bll)  then  gives 


dw 


= (2  sin  1 (2  cos  2 j^(l  + a1  cos  w + a2  cos  2o>  + ...) 

+ (a-|  sin  w + s^n  2u)  + •••) 


2 
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and  the  best  that  can  be  done  with  1^  is  then 

T _ 1 ? (2  sin  a)/2)4~3lJ1  (2  cos  m/2)4~?u2  dt» 

d " c j 2 rm 

0 [(1  + a^  cos  u)  + a2  cos  2u  + ...)  + (a^  sin  w + a2  sin  2w  + ...)  ] 

(BK) 

So  far  no  way  has  been  found  to  avoid  direct  numerical  evaluation  of  this  formula 
except  in  certain  special  cases. 


Obviously  if  the  necessary  number  of  coefficients  is  large  and  the  number 
of  quadrature  points  must  be  large  for  sufficient  accuracy,  the  evaluation  of 
(B14)for  the  many  permutations. required  to  ace  the  minimum  of  F could  become 
a very  tedious  and  costly  operation.  Therefore,  there  Is  a strong  incentive  to 
find  both  more  accurate  (specific)  quadrature  formulas  than  (say) Simpson's  rule, 
and  methods  of  evaluation  at  least  for  A and  c which  do  not  require  further 
quadratures.  The  topic  of  special  quadrature  formulas  is  considered  later 
(sectior,  9b),  but  the  obvious  alternatives  for  c and  A can  be  derived 
immediately  from  (B9)  by  using  the  expansion  in  powers  of  1/c  valid  when 

Id  > i. 


Thus  by  truncating  each  group  in  (B 9)  it  is  a trivial  matter  to  find  the 
triple  product  coefficients  An  In  the  formula 


dz  _ 


(B15) 


by  using  the  polynomial  product  subroutine  mentioned  in  Appendix  A.  Then  it  is 
easy  to  prove  from(B15)  that 


and 


(B16) 
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both  of  which  may  be  useful  computational  formulas  provided  the  convergence  Is 
adequate.  As  will  be  seen  the  convergence  Is  quite  often  poor,  but (B1 6 ) and 
(B17)  do  have  a theoretical  Interest  In  certain  special  cases  quite  apart  from 
their  general  computational  promise. 

The  most  general  case  has  not  been  computed  fully  since  It  was  not  realized 
In  the  early  stages  of  this  study  that  It  would  be  necessary  to  Incorporate  an 
option  providing  for  finite  closure  angles  at  both  ends.  This  feature  greatly 
Increases  the  program  complexity  through  the  poor  convergence  of  the  zero 
factors 


and  the  associated  revised  general  optimization  program  has  not  been  properly 
tested,  although  preliminary  results  have  been  obtained.  This  phase  Is 
described  further  in  section  89  where  it  properly  belongs,  because  most  of  the 
subtleties  necessary  for  the  general  scheme*  were  learned  from  the  p^el imi nary 
study  of  special  cases.  Therefore,  it  Is  both  easier  and  more  logical  to  dis- 
cuss the  special  cases  first  -more  or  less  in  the  order  in  which  they  were 
developed. 

B4.  ELLIPSES 

A natural  starting  point  for  such  optimization  studies  is  the  ellipse 
profile  because  its  characterization  is  so  simple  in  terms  of  the  mapping 
theory.  Furthermore,  not  only  does  it  appeal  to  the  intuition  as  a low  drag 
shape,  but  in  fact  it  was  known  to  have  good  properties  (see  e.g..  Figure  14 
and  section  8 in  this  report). 

In  the  early  stages  it  was  thought  that  perhaps  the  ellipse  could  be 
established  analytically  as  the  optimum  case  with  a continuous  9 (as 
opposed  to  trick  "fairings"),  on  the  grounds  that  the  ellipse  family  is  the 
only  blunt  ended  symmetric  class  (t]  = = it)  which  also  has  every  a 

coefficient  zero  except  a£.  This  conjecture  did  not  turn  out  to  be  correct, 

♦including,  of  coutse,  the  necessity  for  its  introduction  at  all! 
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but  it  was  possible  to  produce  analytic  expressions  for  every  phase  of  the 
ellipse  study  and  so  give  Invaluable  guidance  for  accuracy  requirements  on  the 
later  numerical  work. 

The  essential  formulas  come  from  (69)  by  first  setting  ^ = m2  = 0 

gi  ves 

1 +7  + 7+  - (618) 

as  the  defining  formula  for  all  symmetric  blunt-ended  bodies.  Then  the  ellipse 
family  is  defined  by  further  putting  a3  = a^  = ...  = 0 so 

||  » 1 + ; z = C + c-  ~;  C real 

Choosing  C = 0 to  put  the  origin  at  the  center  and  then  setting  c = e1w 
gives  the  surface  coordinates  as 

xs  = ^ ~ a2^  cos  (1)*  ys  = (1  + a2)  sin  w 


so  that 

a = semi -major  axis  = 1 — a2 
b = semi -mi nor  axis  = 1 + a2 

and  the  range  of  a2  is  -1  <_  a^  1+1.  Clearly  the  likely  (slender)  range  of 
low  drag  shapes  corresponds  to  a2  < 0 and  obvious  limiting  cases  are  the 
streamwise  flat  plate  (a,,  = -1)  and  the  circle  (a2  = 0). 

Now  the  structure  of  ds/dai  from  (B13)  is 

^ ~ + 2a2  cos  2uj  + a2  = ^ja  sin  w + b cos  u> 

from  which  it  Is  clear  that  for  the  range  a2  < 0 (or  b < a)  of  Interest 
here,  integrands  involving  ds/dw  can  be  expressed  In  terms  of  complete  elliptic 
integrals  by  choosing  the  modulus  as 
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thus  giving  the  complementary  modulus  k‘,  where 

k'2  = 1 - k2  = b2/a2 

as  the  fineness  ratio.  Then  the  perimeter/chord  ratio  is 

P ’ c f ■ E 

0 

and  the  fundamental  Ip  Integral  can  be  written 

In  = t(2  " k*>E  " 2k'2^ 

u k^ 

where  E and  K are  the  standard  complete  elliptic  integrals.  Since,  In 
addition. 


these  results  can  be  used  to  calculate  the  behavior  of  the  comparison  factor  F 
of  (B7)  with  great  precision  as  a function  of  the  fineness  ratio  k'  = b/a. 


The  limiting  values  behave  as  expected.  Firstly,  as  k ' -*■  0 (k  1), 

I.  -+■  0 but  K behaves  as  In  4/k'.  The  limiting  value  of  the  last  term  is 

2 

nonetheless  zero  since  it  is  multiplied  by  k'  . Furthermore  E(l)  = 1 so 
the  limiting  value  of  Ig  is  1 as  it  should  be  according  to  the  definition 
(B3)  since  k ' = 0 is  the  flat  plate  case. 


As  k*  1 the  profile  tends  to  circularity  and  I,  -+  tt/8  as  it  should 

* -4 

according  to  (B5).  Some  care  is  needed  with  Ip  as  k -*■  0 owing  to  the  k 
factor,  but  it  turns  out  that 

(2  - k2)E  - 2k'2K  * k4  + 0(k6) 


i 

i 

< 
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giving  the  limiting  value 


l[j(0)  = 3tt 

which  implies  that  a circular  cylinder  in  turbulent  flow  has  a drag  coefficient 
5/6 

(3*)  x flat  plate  of  the  same  chord  - provided,  of  course,  that  separation  is 
ignored 

In  any  event  Fig.Bl  shows  the  ellipse  curve  as  a function  of  fineness  ratio 
(f)  with  the  clearly  defined  minimum  at  about  f = k * =0.22.  Anticipating 
somewhat,  it  can  be  seen  that  this  curve  has  very  nearly  the  lowest  minimum  of 
al 1 those  on  Fig.  B1 . 


B5.  SOME  CUSPED  AIRFOILS 

No  study  like  this  Is  complete  without  the  Joukowski  airfoil.  This  shape 
also  offers  some  hope  of  complete  analytic  evaluation,  but  there  are  others  of 
even  simpler  form.  These  are  all  of  intrinsic  interest  because  the  cusped 
trailing  edge  removes  the  possible  objection  that  any  wholly  blunt  shape  deduced 
as  an  optimum  by  this  Integral  formulation  would  always  be  liable  to  poor 
practical  performance  due  to  separation. 

The  Joukowski  profile  is  a member  of  the  more  general  class  known  as  pole 
airfoils  which  are  mentioned  In  Appendix  A and  discussed  in  more  detail  in  [B2]. 
The  single  pole  airfoil  is,  strictly  speaking,  the  simplest  member  of  this  class 
and  so  will  be  considered  first. 

(a)  A single  pole  airfoil  is  characterized  by  = 0,  = v and 

af  = (’ -tK’ • |d|^'  (B’9> 

d being  real  for  the  symmetric  case.  The  second  factor  can  be  expanded  to  give 
the  an  coefficients  as 
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but  again  (as  a check  case)  a more  direct  evaluation  was  used  for  Ip,  only 
the  area  being  computed  from  the  series.  Thus,  instead  of  expanding  (B19)  the 
modulus  of  1 + 1(s  — d)  can  be  calculated  on  the  surface  and  used  in  (Bll)  to 
give 


2 sin  f 


+ 4(1 



d)  cos_  (»/2 


(1  + d)  - 4d  cos  w/2 


so  that  substituting  in  (B12)  and  using  the  transformation  £ = cos  w/2  yields 


In  - 


64 


/ 


(1  + d)2  - 4dc2 


13/2 


d7  + 4(1 


dk‘ 


54dC 


Since  the  chord  is  given  by 


„_o.l-d-,_/l+d\ 
c _ 2 + — g—  In  (yrra-J 

and  the  thickness  by 

t - 2ymt  ■ /tr-wnr-!^)  sin-'  (42^2) 

computation  was  straightforward,  but  IQ  was  evaluated  by  using  a 101  point 
Simpson  routine  on  looking  at  the  above  formula.  The  results  are  shown  on 
Fig.  B1  and  clearly  exhibit  a minimum  at  about  f = 0.2.  Somewhat  disappoint- 
ingly the  drag  is  evidently  much  greater  than  for  the  ellipse. 


(b)  In  the  pole  format  a Joukowski  profile  taxes  the  form 

af - (i + r^  + ^f^]; 


again  with  d real.  Following  the  same  steps  as  for  the  single  pole  case  leads 
to 

1 

r/l  J.  A\^  AAr^l''  r^Ar 

(B20) 

V AA <-  A / O A 

0 


|4  r 

D c i [4dZ  - 4(2d  - l)cZ]3/ 
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with 


c = 


2d  sin  <i>Q(l  - cos  u>0) 
1 - 2d  cos  to  + 


cos  w. 


= 0 - d + d2)  -Vi  - d2  + d 


2a 


Again  (B20)  was  evaluated  by  a 101  point  Simpson  quadrature,  but  this  time  an 
exact  integration  was  possible  and  was  carried  out.  It  served  to  show  that 
the  results  obtained  by  the  numerical  procedure  were  quite  adequate.  The 
results  are  given  on  Fig.  B1  and  exhibit  a very  similar  behavior  to  the  single 
pole  case  with  a somewhat  higher  drag. 


(c)  The  final  experiment  in  cusped  trailing  edges  was  a single  parameter 
airfoil  which  is.  In  a sense,  the  cusped  analog  of  the  ellipse  since  it  is 
defined  by  all  an  = 0 except  a2  so  that 


(B21) 


On  the  basis  of  the  kind  of  argument  advanced  intuitively  for  the  ellipse, 
namely  that  it  should  be  some  kind  of  optimum  since  all  the  higher  order  modu- 
lations due  to  an  are  absent,  it  might  be  expected  that  (B21)  would  be 
superior  to  both  the  Joukowski  and  single  pole  profiles  since  they  both  have 
an  infinite  sequence  of  an  coefficients. 


However,  some  doubt  can  be  cast  on  this  Idea  immediately  because  the 
modulating  g(c)  factor  in  (B21)  is 


g(c) 


c2  + c 


+ a. 


where 


cl*  H 


If  1 1/4  both  zeros  are  real  and  since  no  zero  can  lie  outside  the  unit 
circle  the  overall  permissible  range  of  a2  Is  0 £ a2  <.  1.  When  the  roots 
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become  complex  the  profile  starts  to  bulge  and  unlike  all  the  previous  cases 
the  limiting  form  at  a2  = 1 has  tw0  symmetrically  disposed  cusps  just  down- 
stream of  the  leading  edge!  Consequently,  it  is  clear  that  (even  without 
separation)  this  shape  will  have  very  poor  drag  characteristics  for  a2 
appreciably  greater  than  1/4, 


Using 


Forwarned  by  this  argument  the  numerical  result  comes  as  no  surprise. 

c - 4(1  -}a2) 
t _ 3/T  , 

t - ~ra  2 


and  the  original  form  of  Ip  (i.e.,  (B14))  with  a Simpson  quadrature  gives 
the  curve  shown  on  Fig.  B1 . Evidently  the  characteristics  are  much  worse  than 
the  previous  cusped  cases,  which  serves  to  show  very  clearly  the  naivete  of 
the  conjectures  made  before  about  optima  coming  from  having  as  many  an 
coefficients  as  possible  vanish.  Further  very  surprising  evidence  concerning 
this  situation  is  provided  by  the  next  section. 


B6.  INITIAL  OPTIMIZATION  FOR  BLUNT  BODIES 


The  evidence  of  section  B5  lent  weight  to  the  idea  that  the  really  good 
shapes  would  have  fore  and  aft  symmetry.  Furthermore,  there  was  at  least  a 
suspicion  that  the  cusp  was  an  undesirable  feature  quite  apart  from  its 
unsymmetrlcal  role  in  the  airfoils  studied.  Therefore,  it  was  natural  to  con- 
sider the  problem  of  minimizing  F only  for  blunt  ended  shapes.  Under  these 
circumstances  ■ n,  t2  = n so  that  = 0,  u2  = 0 and  the  fundamental 
formulas  become  simplified  mainly  because  then 


so  that  (B16)  and  (BI7)  could  be  used  directly  without  undue  fears  about 
convergence.  In  fact  (B7),  (B8),  (B14),  (B16)  and  (B17)  show  that  F is  a 
function  only  of  the  an  numbers  for  blunt  ended  shapes,  and  a definite 
cut-off  for  the  an  sequence  consistent  with  some  optimistic  belief  about 
convergence  is  obviously  a necessity  for  any  numerical  work. 


103 


On  this  basis  a computer  program  was  written  to  effect  the  minimization, 
and  it  can  be  described  very  simply  as  follows: 

(a)  It  was  decided  that  rather  than  try  sophisticated  gradient,  steepest 

descent,  etc.  methods,  a very  simple  "boxing"  procedure  to  progressively  narrow 

iri  on  the  minimum  of  F(a|,  a2>  ...  an)  would  be  adequate,  mainly  because  of 

the  conviction  that  it  only  needed  to  be  done  once  in  order  to  isolate  the 

optimum.  This  very  simple  program  works  by  the  following  steps  which  depend 

only  on  a subroutine  to  evaluate  F qiven  a set  of  a„  values. 

n 

Step  1.  Input  an  initial  guess  for  the  an  values  and  a starting 
step  size  h. 

Step  2.  Increase  a-|  to  a-j  + h.  If  F has  decreased,  move  on  to 
a2.  If  F has  increased,  try  a-j  - h.  If  this  also  gives 
an  increase,  leave  a^  alone  and  move  on  to  a^.  Continue 
this  process  to  an> 

Step  3.  After  one  pass  as  in  step  2,  repeat  the  cycle  to  see  if  the 
changed  values  of  a^,  a^»  ...  have  changed  the  decision 
regarding  a^ . Continue  this  checking  until  acomplete  pass 
does  not  change  F. 

Step  4.  Decrease  h by  a factor  of  10  and  start  again  at  step  2. 

Step  5.  When  h has  been  reduced  to  desired  accuracy,  stop. 

It  can  be  seen  that  each  completed  cycle  (end  of  step  3)  means  that  to  within 
the  current  value  of  h,  it  is  not  possible  for  the  program  to  detect  a 
further  decrease  in  F,  which  implies  that  the  local  minimum  is  contained 
within  a hypercube  (box)  of  volume  hn.  The  division  of  h by  10  each  time 
is  a particularly  convenient  mode  of  operation  since  then  each  completed  cycle 
defines  another  decimal  place. 

Various  runs  were  tried  to  test  this  program  using  special  test  functions 
of  several  variables  and  the  direct  results  already  reported  for  F depending 
on  only  one  variable.  It  was  very  soon  appreciated  that  high  accuracy  for  the 
quadratures  was  necessary  if  a large  number  of  terms  were  to  be  Isolated, 
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since  the  whole  process  hinges  on  the  detection  of  a genuine  change  of  F due 
to  a change  of  an  as  opposed  to  variations  caused  by  numerical  error.  This 
topic  will  be  discussed  later  (Section  D9);  it  is  sufficient  to  record  here 
that  50  - 101  quadrature  points  were  adequate  for  20  - 30  an  coefficients 
with  h as  small  as  0.00001.  For  this  particular  version  of  the  program 
further  description  is  not  necessary  because  of  the  results  that  it  led  to 
immediately. 


Of  course,  in  the  case  of  the  ellipse  it  was  not  necessary  to  worry  about 
accuracy  of  quadrature  formulas  since  all  the  necessary  calculations  were 
exact.  Under  these  circumstances  it  was  a trivial  exercise  to  determine  the 
optimum  ellipse  by  the  minimization  program  and  the  result  when  cycled  down 
to  h = 0.0000001  was 

a2  = -0.6375691,  f = t/c  = 0.2213225,  F = 7.7122530 
This  is  shown  on  Fig.  B2  compared  with  some  other  "good"  profiles. 


( b)  The  first  two  serious  attempts  to  find  general  optima  were  with  10 
and  20  an  coefficients.  Since  the  result  obtained  with  the  20  term  case  was 
a repeat  of  the  10  term  with  greater  accuracy,  we  shall  not  discuss  the  10 
term  at  all.  In  fact,  both  runs  were  first  carried  out  without  the  5/6  power 
on  Ip,  but  this  only  changes  the  details  not  the  nature  of  the  answers. 

First  It  was  observed  that  as  h was  decreased  through  each  cycle  all  the 
odd  coefficients  decreased  rapidly  which  implies  very  clearly  that  fore  and  aft 
symmetry  Is  a necessary  requirement.  After  5 cycles  (h  = 0.00001)  the  follow- 
ing even  coefficients  were  obtained: 

a2  = -0.71513 
a4  = -0.09447 
a6  = -0.04134 
ag  = -0.02384 
cl i q - -0.01567 
a12  = -0.01111 
a-|4  = -0.00825 
a16  = -0.00628 

a18  = "°*00480 
a2n  = -0.00349 
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It  was  noted  that  this  body  did  have  a slightly  lower  F than  the  best  ellipse 
by  about  G.15%,  but  by  far  the  most  striking  feature  of  the  coefficients  was 
their  poor  convergence.  On  plotting  the  shape  it  was  noted  that  the  "bluntness" 
was  confined  to  an  exceedingly  local  region  near  the  ends,  with  the  overall 
aspect  having  every  appearance  of  a finite  angle  wedge  ended  body.  By  a very 
crude  protractor  estimate,  a value  of  about  50°  was  obtained  from  the  plot 
leading  to  t-j/n  % 0.28  or  ^ = y2  ^ 0.72.  In  view  of  the  leading  coefficient 
in  the  above  table,  the  obvious  requirement  for  fore  and  aft  symmetry  and  the 
poor  convergence,  it  was  only  natural  to  compare  the  table  values  with  the 
expansion  of 


for 


l 

7 


0.72.  The  values  of  these  coefficients  are 

a2  = -0.72 
a4  = -0.1008 
3g  = -0.0430 
ag  = -0.0245 
a1Q  * -0.0161 
a12  = -0.0115 
a14  --  -0.0087 
a1g  = -0.0068 
a^  = -0.0055 

**20  = 


On  comparing  the  two  sets  of  values  and  remembering  that  the  "minimization" 
sequence  was  truncated  and  that  the  estimate  for  t was  very  crude,  it  was 
hardly  possible  to  avoid  the  Inference  that  the  minimization  procedure  was 
leading  to  a limiting  form  of  solution  described  by 


dz 

3? 


('-W* 


u ^ 0.72 


Inasmuch  as  such  bodies  would  have  finite  trailing  edge  angles  and  character- 
istic rates  of  convergence  as  shown  above,  this  discovery  was  singularly 
unwelcome  as  well  as  surprising.  It  meant,  at  least,  that  much  mere  serious 
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thought  would  have  to  be  given  to  the  structure  of  subroutines  designed  to 
evaluate  F for  use  In  the  minimization  program.  However*  this  Is  once  again 
best  discussed  in  Section  B9,  since  It  is  a topic  relevant  to  the  most  general 
program  and  the  important  details  can.  In  any  case,  be  best  appreciated  by 
examining  the  simpler  special  solutions  with  finite  closure  angles  first. 

B7.  CANONICAL  WEDGE  END  PROFILE 

If  the  ellipse  can  be  considered  as  the  canonical  blunt  end  shape  on  the 
grounds  that  it  has  only  one  an  coefficient,  then  the  class  suggested  in 
Section  B6  as  the  tentative  limit  of  the  minimization  procedure,  namely 

(l  _l_ju  , 0 < M < 1 (B22) 

c 

can  certainly  be  regarded  as  the  canonical  wedge  end  profile  since  it  has  none. 
Once  again  the  temptation  exists  to  jump  to  the  conclusion  that  because  there 
are  no  an  coefficients  this  profile  must  be  the  optimum;  but  having  been 
wrong  twice  before  (section  B4  and  B5c)  it  is  possible  to  defer  judgment  on 
this  issue. 

The  consequences  of  (B22)  are  quite  interesting.  Since  all  the  an 
coefficients  are  zero  the  expansion  of  (B22)  gives  the  Ap  coefficients 
immediately  as 


a = h.(1  Id. 

n n! 


so  that 


n+1 


= k -Hr 


and  therefore  the  expansion  is  uniformly  convergent  [B3]  when  |c|  >_  1 if 

u > 0 or  n > t 

which  covers  all  cases  of  interest  here.  Hence  (B16)  and  (B17)  are  apparently 
trivial  calculations  to  give  area  and  chord.  Furthermore,  from  (BIO)  and  (Bll) 
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and  so 


(2  sin  w)1J; 


Qs  - (2  sin  u)T/v 

which  are  all  unusually  simple. 


However,  the  most  Important  feature  is  that  the  drag  integral  also  turns 
out  to  be  an  analytic  result.  This  can  be  seen  from  the  combination  of  the 
above  results  because 


7T  TT 

1 = J Q4  dw  = j (2  sin  u)3e+1  dw 
0 0 

where  we  have  used  e for  x/ir  as  a convenience  so 

p = 1 — e 


Using  <j>  - v/2  — u this  can  be  written  as 


,3e+l 


7T/2 


3e+1 


•I  = -2 

and  then  by  using  the  expansion*  (-tt/2  <_  <p  <_  tt/2) 


f (cos  <i> )''c' ' ’ d<f> 

in 


cos”*  = —V , Ti  ♦ -fy  COS  24>  + r 

2V"'  r[()/2)v  + l]2  L?  v + z Tv  + 27Tv 


IV  cos 


t ...j 


it  follows  that 


I = r(3e  + 2Jtt 

f[(3e/2)  + (3/2)]* 


ID  - I/c 


another  remarkably  simple  result.  It  can  be  seen  that  this  has  the  correct 
behavior  in  the  limits  since 

2 

a.  e = 0 gives  I = r(2)Tr/r(3/2)  = 4 which  is  correct  for  a mapped 

flat  plate  which  always  has  a chord  of  4 giving  Ip  = 1. 

*This  appears  In  [B3],  p.  263,  No.  40.  It  is  also  uniformly  convergent  when 
e < 1 and  has  the  same  asymptotic  rate  of  convergence  as  the  Ap  expansion. 
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p 

b.  fc  = 1 gives  I = r(5)n/r(3)  *=  6tt  so  1^  = 3n  since  c = 2 

for  the  circle  which  Is  the  correct  result  according  to  section  B4. 

By  using  these  formulas  the  analytic  behavior  as  recorded  on  Fig.  B1  was  easy  to 
calculate. 

However,  this  calculation  served  to  Indicate  a rather  obvious  weakness  in 
the  expanded  forms  used  for  the  chord  and  area.  In  order  to  determine  the  chord 
accurately  a lot  of  terms  had  to  be  used  since  the  convergence  of  the  An 
sequence  is  poor.  To  some  extent  this  was  also  true  of  the  area  but  much  less 
severe.  Thus,  at  a f=  t/c  ratio  of  0.215  it  was  necessary  to  use  200  terms 
In  (B1 6 ) ! 

The  subroutine  which  generated  the  various  formulas  was  also  cycled  through 
the  optimization  program  described  before  (section  B6a)  down  to  a step  size  of 
h = 0.0000001,  and  gave  the  optimum  case  of  this  family  as 

e = 0.3397098,  f = t/c  = 0.2151623,  F = 7.7007989 

In  this  application  - as  in  the  straightforward  parametric  variation  to  define 
the  curve  of  Fig.  Bl  - the  slow  convergence  was  not  important  since  only  one 
parameter  (u  or  actually  e)  was  being  varied,  but  the  implication  for 
general  studies  with  many  parameters  were  not  regarded  as  very  encouraging. 

This  will  be  discussed  further  in  section  B9.  For  the  moment  we  merely  note 
that  the  optimum  member  of  this  family  (above)  has  a wedge  half  angle  of  about 

0.3397098  x 90  * 30.57° 

and  is  shown  on  Fig.  B2.  According  to  section  B6a,  the  value  of  F is  certainly 
lower  than  that  of  the  opt^ium  ellipse,  but  the  gain  Is  trivial  since  the  CDa 
ratio  is  only  1.0015  (see  Fig,  B2). 

B8.  CIRCULAR  ARC 

Because  of  the  evident  importance  and  different  nature  of  wedge-end 
profiles  it  was  decided  that  at  least  one  other  case  should  be  Investigated. 

It  is  known  that  the  profiles  generated  by  segments  of  circles  can  be  analyzed 
exactly  and  in  fact  the  basic  solution  is  given  by  Milne-Thomson  [B4].  However, 
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one  of  the  reasons  for  considering  this  case  was  to  compare  the  known  theory 
with  the  mapping  format  used  here  - In  order  to  get  a better  understanding  of 
the  convergence  situation.  This  was  considered  necessary  because  section  B? 
shows  how  poorly  the  zero  factors  converge  when  expanded  naively,  but  this  was 
no  guarantee  that  a modulating  term  containing  an  coefficients  would  not  make 
the  situation  worse  still.  Obviously  such  Information  and  test  cases  would  be 
important  In  designing  an  overall  subroutine  suitable  for  the  most  general 
minimization  applications. 

Therefore,  some  detail  had  to  be  added  to  Milne-Thomson's  exposition  and 
this  turned  out  to  be  very  fruitful  since  an  exact  formula  (no  quadratures)  was 
found  for  the  drag  integral.  This  case  then  became  the  most  important  test 
problem  since  it  had  all  the  features  desired  for  the  most  general  program  and 
all  the  relevant  formulas  exact. 


Hence  some  explanation  cannot  be  avoided  and,  whilst  circular  arc  profiles 
are  a natural  aspect  of  Karman-Trefftz  mappings,  the  coaxial  circle  method 
of  Milne-Thomson  is  preferrable  since  it  (eventually)  suggests  the  forms  which 
permit  exact  evaluation  of  the  l)^  integral.  Briefly,  the  classical  procedure 
is  to  start  from  bipolar  coordinates  defined  in  the  z = x + iy  plane  as 
sketched,  whence 


so  that  profiles  of  constant  x are  circular  arcs.  The  mapping 

z = a coth  ^ v (B23) 

then  maps  the  lines  -»  <_  u < +»,  0 >_  \ -2n  into  a connected  sequence  of 

shapes  consisting  of  the  real  axis  and  a piece  of  circular  arc  which  can  either 
be  concave  or  convex  from  a chosen  side. 
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(a)  There  are  many  ways  to  get  fluid  motions  Into  such  a system  (see  e.g., 
Mllne-Thomson),  but  one  of  the  most  fruitful  Is  to  consider  various  W -*■  v 
mappings  where 

W = <f>  + 1#  = complex  potential 

In  particular  the  relation 

“ ■ r “«*  (i) 

is  Interesting  because  separating  real  and  Imaginary  parts  gives 

* - («*) 


the  first  of  which  vanishes  for 

x = 0 , +_  5p-  , +_  mu,  ... 

Leaving  aside  some  slightly  subtle  details  about  branches.  It  can  be  shown  that 
if  0 m 2 then  the  cases  X ~ 0,  X = - m/2  give  flows  having  the  real 
axis  and  one  of  the  above  sequence  of  circular  arcs  as  streamlines. 


Although  it  is  not  possible  to  eliminate  v directly  to  give  a W(z) 
relation  in  general,  the  two  limiting  cases  of  interest  to  this  work  can  be 
done,  viz. 

when  m = 2,  w = z (infinite  plate  along  y = 0) 

p 

when  m = 1 , w=z+a/z  (semicircle) 
and  evidently  the  range  of  interest  is  1 _<  m 2. 


(b)  The  complex  velocity  is,  however,  expressible  in  terms  of  v for  all 
cases  since 


q = u 


iv 


dW  = dW  dv 

37  ~ dv  3z 


4 sinh2  v/2 

T — 2 
m sinh  v/m 


so  that  the  surface  velocity  magnitude  can  be  written 


« _ 4 / cosh  u - 1 \ 

% ^ \cosVW-  1 ) ’ 

- 4 / cosh  \i  - cos  rmr/2  \ 

~7  \ cosh  £u/m  + 1 / 

which  are  reasonable  enough  formulas  in  terms  of  the  parameter  u.  However, 
to  relate  these  to  the  unit  circle  plane  contained  in  ? = z + in  it  is  neces- 
sary to  find  a mapping  which  maps  the  unit  semicircle  into  one  of  the  segments 
and  has  the  same  distant  regions.  Evidently  a mapping  of  the  same  kind  as 
(B23)  must  suffice,  and  by  considering,  therefore 

t,  = c coth 

it  is  easy  to  show  that  the  choice  n = 2m,  c = 1 Is  required  giving 
' = coth  £ or  e'/m 

Comparing  this  with  (B24)  and  noting  that  the  distant  fields  match  only  when 
m = 2a  gives  the  relation 

z + m _ i r.  + 1 ; m = 2 — e = 1 + u (B27) 

z — m \ c - 1 ' 


2 2 
x > a 


-a  < x < +a 


(B26) 


which  is  the  usual  Karman-Trefftz  form. 

By  means  of  (B26)  and  (B27)  the  surface  velocity  can  be  expressed  in  terms 
of  the  usual  unit  circle  angular  variable  since  we  get 


u/m  _ 


= cot 


and  then 


•j 

n - * to  /_2m  b)  . » 2m  a)  A o Tf  __m  to  _i_m  to  \ 

Qs  - — ig-  (2  sin  w)  (cos  -g  + S1n  y + ^ cos  e 2 cos  S’  s'n  2 I 

m ’ 

which  is  fine  for  specific  computation  of  Q$  on  the  profile,  but  not  of  much 
general  use. 
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Much  more  significant  Is  that  by  using  <*>  as  the  Independent  variable  an 
exact  form  for  the  basic  drag  Integral  can  be  found.  Thus,  from  (B25)  with 
X = -mir/2  and  m = 2a 


giving 


♦ 


2 sinh  2p/m 
cosh  2p/m  + 1 


2 tanh  £ 
m 


2+T 

r=t  * 


(-2  <_  <()  £ 2 on  the  surface) 


whilst 


i / /q4| /q3^ 

0 -2  -2 


and  Q can  be  written  (with  0 = cos  ett/2) 


(B28) 


Q = [(2  + 4>)m  + (2  - <f>)m  + 20(4  - *2)m/2]  (4  - *2)c/2  (B29) 

4m  L J 


Substituting  (B29)  Into  (B28)  and  accounting  for  many  symmetries  In  groups  of 
terms  gives,  finally,* 


1 = fr^7  + T*  + 6sr<6  “ + E)  + 30  + 4ti2)r(5  -f)r(3  + 4) 

m r(8)  L c £ * * 

+ 26(3  + 202)r(4)r(4)J  (B30) 

The  other  quantities  needed  for  direct  evaluation  are  the  chord  and  area.  These 
are  given  very  simply  by 

c = m (B31 ) 

and 


*It  is,  again,  easy  to  show  that  (B30)  gives  the  same  correct  limiting  values  of 
Iq  in  the  flat  plate  and  semicircle  cases  as  before,  but  there  Is  no  point  in 
giving  the  details  here. 
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where  .)  is  half  the  closure  angle  (e  = T/2).  Note  that  for  small  values  of 
0,  A -*•  0 of  0(e),  in  fact 

A ^ m2  [f  0 + ^5-  e3  + 0(e5)] 

Putting  (B30),  (B31)  and  (B32)  together,  the  circular  arc  curve  of  Fig.Bl  was 
produced  by  completely  analytic  steps.  These  formulas  were  also  run  through 
the  minimization  program  down  to  h = 0.0000001  with  the. resulting  optimum 
given  as: 

e = 0.2693423  (t  * 48. 5°),  f = 0.2147540,  F = 7.7239339 

As  can  be  seen  from  Fig.Bl  the  drag  characteristics  are  (slightly)  less  desir- 
able than  either  the  ellipse  or  canonical  wedge  end  profiles.  Fig.B2  shows  the 
optimum  shape  compared  with  the  optima  of  these  other  two  families  and  also 
indicates  that  the  drag  ratio  of  the  best  circular  arc  is  only  1.0030  x that 
of  the  best  overall  wedge  end  case. 


(c)  However,  as  pointed out  before  the  greatest  value  of  the  circular  arc 
is  as  a test  case,  because  it  has  all  the  desired  features  of  the  general 
method.  It  remains  therefore  to  cast  the  unit  circle  mapping  (B27)  into  the 
format  of  this  work.  Differentiating  gives 


dz  _ 4m2  (C  + l)M(c  -1)M 

^ lu  + i)m-u-i)mr 


u = l--e  = m — 1 


which  can  be  reorganized  to 


2m 

SCO  + 1/C)m  - (1  - l/cf] 


showing  that  the  modulating  g(^)  part  is 


g(  c ; a ~^K 

c 


[o  +i)"-n  -^)m]"2 


p 

Thts  can  be  expanded  In  powers  of  1/s  since 
1 *m 


where 


so  that 


O +i)l-(l  = (l  +a2  ^ 


+ ^/l  T + • • • 


47 


a,  = >ii Hg-.U  , ^ = hLk.-J .),(,P,c72)(H  -3) 


i(0  ■ (i 


a2  A a4 


g(0  - (i  +-f  + -r  + ...)  = 1 - 

c c 


-2 

•) 


5! 


2a2  (2a4  + 3a2) 

~7 


These  coefficients  can  be  computed  to  any  reasonable  order  by  using  existing 
subroutines  which  operate  on  series  and,  as  expected,  g(0  Itself  does  not 
converge  impressively.  For  Instance,  with  v * 0.9  (which  Is  certainly  part 
of  a reasonable  test  range  ) the  40th  coefficient  of  (1  + l/c)p  Is  0.00008737 
and  the  40th  coefficient  of  g(0  is  0.00004556,  30  the  rate  of  convergence  is 
about  the  same  as  for  the  singular  factors.  The  combined  40th  term  is  0.00034836 
showing,  as  anticipated,  that  product  of  the  poorly  convergent  singular  factors 
with  the  g(c)  sequence  only  makes  things  worse.  All  this  information  is 
germane  to  the  structure  of  the  general  minimization  program.  This  topic  Is 
discussed  next. 


89.  GENERAL  MINIMIZATION  PROGRAM 

(a)  The  numerical  work  described  in  sections  B4  - B9  (Fig. B4)  succeeded 
in  establishing  that  the  canonical  wedge  end  optimum  (section  B7)  was  clearly 
the  best  of  all  the  families  and  optima  considered.  However,  as  remarked  before 
on  a number  of  occasions  (e.g.,  section  B5c)  we  were  loath  to  jump  to  the 
conclusion  that  this  was  the  optimun  just  because  of  Its  extreme  simplicity. 

Not  only  were  such  hasty  arguments  disproved  by  further  calculations  before 
(e.g.,  section  B5c),  but  also  it  should  be  noted  that  the  circular  arc,  which 


115 


was  Introduced  merely  as  a further  check  case  and  which  has  a complicated  and 
poorly  convergent  structure  in  terms  of  the  an  coefficients,  nevertheless 
gives  an  optimum  very  close  to  the  canonical  shape.  Moreover  another  factor 
had  entered  the  picture  at  this  point.  The  complex  calculus  of  variations 
theory  had  been  completed  - admittedly  on  a very  tentative  basis.  However, 
neither  the  ellipse  (special  case  for  t = tt)  nor  the  canonical  wedge  end 
profile  satisfied  the  resulting  differential  equations  supposedly  defining  the 
real  extremum*.  Although  not  conclusive,  this  fact  together  with  the  above 
arguments  was  sufficient  to  indicate  that  further  numerical  minimization 
studies  were  needed,  and  that  they  must  encompass  the  complete  system  as  out- 
lined originally  in  section  B3. 

As  has  been  noted  before  the  convergence  for  the  series  In  the  general 
case  is  not  adequate  for  them  to  be  used  to  determine  the  area  and  chord. 
General  operation  of  the  numerical  scheme  described  In  section  B6a  requires 
a large  number  of  passes  for  (say)  the  determination  of  20  - 40  coefficients, 
and  procedures  adequate  for  only  one  parameter  have  to  be  reconsidered  owing 
to  the  computing  costs  — even  though  this  particular  calculation  only  needs  to 
be  done  once  in  principle.  In  this  context  It  is  also  worth  noting  that  if 
this  is  true  in  two  dimensions.  It  is  obviously  even  more  Important  for  the 
axisymmetric  case. 

(b)  At  first  it  was  hoped  that  the  convergence  difficulty  could  be  circum- 
navigated by  replacing  the  series  for  chord  and  area  by  quadratures.  A new 
program  was  written  on  this  basis  using  various  different  numerical  methods: 
Simpson's  rule 
Cubic  spline  fits 
Fast  Fourier  transform  fits 

to  carry  out  the  quadratures.  Eyen  using  201  points  not  one  of  these  was 
remotely  accurate  enough  when  compared  with  the  circular  arc  exact  solutions 
in  the  region  of  the  various  optima. 


♦No  solutions  for  these  equations  have  been  found  ab  Initio  partly  because  of 
the  speculative  nature  of  their  derivation.  This  phase  of  the  work  Is  simply 
not  developed  enough  at  the  time  of  writing. 
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The  reason  for  this  failure  of  the  conventional  methods  to  give  the  kind 
of  precision  necessary  for  the  minimization  process  to  detect  the  Influence 
of  higher-order  small  coefficients  Is  not  hard  to  understand.  All  of  such 
conventional  methods  depend  very  strongly  on  the  fundamental  axiom  of  classical 
numerical  analysis  — that  the  data  corresponds  to  a sufficiently  smooth  (differ- 
entiable) function.  However,  the  functions  to  be  integrated  here  are  of  the 
form  (e.g,,  see  (814)). 

TT 

I = £ (2  sin  |)P  (2  cos  |)q  f(u)  dw  (B33) 

where  f(w)  is  smooth  in  the  desired  sense,  but 

0 i p.q  1 1 

Thus  the  singular  factors  cause  a very  rapid  rise  at  the  ends  of  the  range  and 
in  fact  the  Integrand  of  (633)  behaves  as  J3  for  small  u and  so  has  an 
infinite  derivative  as  w -*■  0. 

One  way  to  get  around  this  difficulty  is  the  old  (but  very  Important) 
artifice  of  separating  the  singular  terms  and  doing  them  analytically.  A sub- 
routine to  affect  this  has  been  written  and  tested  for  Integrals  of  the  type 
shown  In  (B33).  Its  efficacy  Is  staggering  as  can  be  seen  from  the  results 
tabulated  below  on  a simple  test  case  where  the  complete  answer  Is  known,  for 
Instance  f(w)  = 1. 

(c)  The  general  program  using  this  artifice  for  nunerlcal  quadrature  has 
not  been  written  but  would  not  be  difficult  since  the  subroutines  required  and 
the  basic  operating  process  of  the  minimization  part  are  all  available. 

However,  this  program  would  not  just  be  a repeat  of  that  of  section  B6a  with 
more  accurate  quadrature,  since  another  very  Important  feature  for  time  and 
cost  saving  Is  that  the  old  program  computed  the  whole  Integrand  each  time, 
whereas  It  Is  obviously  only  necessary  to  update  one  term  in  (e.g.)  the  Fourier 
expansions  each  time  an  an  coefficient  is  changed. 
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Table  Bl. 


Error  of  I for  Various  Methods  When  f( 


Simpson's 

Rule 


-0.390048958 

-0.224950593 

0.127376081 

0.070438403 

0.037729564 

0.019315228 

0.009221491 

0.003888859 

0.001220631 

0.000000540 


-0.396812512 

-0.232855499 

-0.134802905 

-0.076724170 

-0.042441076 

-0.022617007 

-0.011301525 

-0.004997188 

-0.001638814 

0.000002570 


1 

Fast 

Fourier 

Transform 

Separation 

of 

Singularity 
+ Simpson's 
Rule 

-0.45374470 

0.000000076 

-0.267544643 

0.000000236 

-0.155229725 

0.000000398 

-0.088178328 

0.000000521 

-0.048640337 

0.000000596 

-0.025747832 

0.000000620 

-0.012755969 

0.000000602 

-0.005608282 

0.000000546 

-0.001845777 

0.000000479 

0.0 

0.000000041 

Since  this  would  not  be  very  difficult  or  costly  it  seems  that  it  should 
be  done  if  only  to  give,  once  and  for  all,  the  final  optimum  - and  so  define 
in  a watertight  manner  the  absolute  lower  limit  on  Cq^.  Not  only  will  such 
studies  not  need  to  be  repeated,  but  engineering  compromises  could  then  be 
organized  on  a firm  objective  basis. 

B10.  CONCLUSIONS 

The  studies  carried  out  as  described  in  sections  B4  ■ B9  lead  to  the 
following  straightforward  conclusions: 


A true  optimum  exists  and  can  be  found  by  a revised,  more  general 
version  of  the  minimization  program. 

There  Is  very  little  to  choose  between  the  class  optima  of  ellipses, 
circular  arcs  and  canonical  wedge  ends. 

However  the  canonical  wedge  end  is  the  best  of  these  and  It  has 
finite  (wedge)  angles  at  the  ends. 
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d.  It  would  not  be  difficult  to  construct  and  run  the  general  program 
(see  a)  to  decide  whether  or  not  c is  the  absolute  optimum. 

e.  Since  the  axisymmetric  drag  integral  is  similar,  it  should  be 
possible  to  repeat  this  procedure  using  the  theory  of  Appendix  A 
and  so  define  an  axisymmetric  true  optimum. 
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